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THERMAL- STRESS ANALYSIS FOR A WOOD COMPOSITE BLADE 


Introduction 

Since 1977, NASA Lewis Research Center has been pursuing the development 
of low cost rotor blade technology for large horizontal axis wind turbines* 

This effort is a part of the Federal Wind Energy Program under the direction of 
the u- S. Department of Energy. Laminated wood manufactured by bonding together 
1/1 0" to 1/8" thick sheets, or plys, with epoxy is a particularly attractive 
candidate material for a rotor blade, since the raw material has high specific 
strength, high specific stiffness and low cost* This led to the manufacturing 
of eight laminated wood blades in 1980. They were installed on Mod-OA wind 
turbines [1] at three different locations; Kahuku Point, Hawaii, Culebra, 

Puerto Rico, and Block island, Rhode Island. Mod-OA turbine has 125* diameter 
rotor and generates 200 kW of power when in operation. Their performance was 
closely monitored. In 1981, a crack was found in Blade No. 1012 located in 
Block Island. The crack occurred in the leading edge extending along the 
entire blade. The cause of the crack was unknown, but thermal-stress induced 
by solar heating was suspected to be one of the main reasons. Therefore, the 
investigation has been pursuaded along this line. 

The analysis of thermal-stress of the blade induced by solar insolation 
consists of two phases. The first phase is to find the temperature distri- 
bution throughout the blade, a heat conduction problem. The second phase 
is to determine the thermal-stress distribution of the blade caused by the 
temperature distribution found in the first phase, a thermal-stress analysis 
problem * Since wood is an ortho tropic material in the senses of heat 
conduction and stress-strain relationships, these characteristics must be 
included in the analyses of both phases. 


Phase I: Heat Conduction Ana lysis 


A. The Blade and the Assumed Environment 

Blade No. 1012 is 700 inches long with maximum chord width of 62*4 inches, 
see Fig. 1 for its overall dimensions. Airfoil standard section NACA 230xx 
series [2] , is used* The thickness to chord ratio of the blade is varied from 
0.3175 at station 150 to 0.075 at station 732* 

In the determination of temperature distribution throughout the blade, it 
is reasonable to assume that a two-dimensional heat conduction analysis of a 
typical blade section is adequate, because the blade is basically a slender 
body. In addition, the thermal conductivity of wood in the longitudinal 
direction of the blade, which is parallel to the grain, is two to three times 
higher than those in transverse directions. As a result, the temperature 
gradient in the longitudinal direction is far smaller than those in the 
transverse directions. This again indicates that a two-dimensional analysis 
is adequate. 

The section at station 126, see Fig* 2, is arbitrarily selected as a 
typical section for analysis. The thickness to chord ratio at this station 

i 

is 0.304, which is close to those of the sections in the main portion of the 
blade. Thus, it is a good representation. In the nose portion, called the 
D-spar, 3" thick laminated rotary veneers of Douglas fir is used. It is 
covered with a layer of 1/8" thick birch plywood. The tail portion is a 
sandwich construction to provide stiffness and light weight* It consists of 
a 3/4" thick honeycomb paper core between two layers of 1/8" thick birch 
plywood. Pieces of sawn Douglas fir are used for stringers. 

Two parked positions of the blade are considered; one keeps the section in 
horizontal position, the other, vertical. Technically speaking, the difference 


between these two positions is that the heat flux input zones are different as 
shown in Fig. 3 and 4. The heat flux input zone is defined as the portion of 
boundary surface of the blade exposed to the sun. The solar insolation is 
estimated at 363 BTU/f t 2 -*hr £ 3 1 • Th e absorbtivity of the blade surface is 
estimated at 0*9, which is not only on the safe side for the purpose of stress 
analysis but also has a good reason: after a long period of exposure, the 

blade surface is usually dirty and full of bug remains, which often drastically 
increase the surface absorbtivity* Thus, the heat flux input, q r , for the heat 
flux input zone is 

q r = 363 X 0.9 = 326 BTU/ ft 2_ hr . 

The ambient temperature, T , is assumed to be 90° F, a typical temperature for 

CO 

a hot summer morning in Block Island. Further, the air is assumed to be 
still, and the convective heat transfer coefficient, h, for the blade surface 
is assumed to be 5 BTU/f t 2 _hr.°F ■ 

In the desert ption of wood properties, three mutually perpendicular axes, 
longitudinal, racial and tangential, need to be used. The longitudinal axis, 

L, is parallel to the grain; the radial axis, R, is normal to the growth rings; 
and the tangential axis, T, is perpendicular to the grain but tangent to the 
growth rings. The thermal conductivities of the wood in L, T and R directions 
[4] are listed in Table 1 f note that a small difference of 10% is assumed 
between k T and to account for possible physical differential between 
these two directions. Since rotary veneer is assumed to be used in blade 
construction, the axes h, T, and R correspond to principal material axes £ , £ 
and r\ respectively. The material axes will be used in next section. 

The thermal conductivity of honeycomb paper core may be estimated by 
using the average value for paper and air, which is 0*013 B TCJ/hr.ft. °F * 




Table 1 . Thermal Conductivity of Wood 


Wood 

k L 

krp 

k R 

Douglas Fir 

0.172 BTU/ hr . ftt o F 

0 . 068 BTCJ/ hr • f t » 0 F 

0.075 BTU/j^. £ t. op 

Birch 

0.211 BTU/ hr . ft ,o F 

0*084 BTU/hj* ^ ^ ^ ^ 0 p 

0.093 BTU/ hr>fti o F 


Bp Governing partial Differential Equation 

Because of the conf igurational complexity of the blade section and the 
orthotropic property of heat conduction of the wood composite, the finite 
element method with variational approach is chosen as the method of solution* 
This method may also be called as a numerical method in solving partial 
differential equations. For an element of an ortho tropic material, the 
controlling partial differential equations may be derived in the following 
outlines. 

Heat flux input - Heat flux output = Heat retention in body ( t ) 

Referring to Figure 5, let *f f be the heat flux and 'k', the thermal 
conductivity, their directions along axes y, 5 and r\ are indicated by 
subscripts* ’T* represents the temperature, ’Cp 1 , specific heat and *p 1 , 
the density of wood. Thus, equation (1) may be expressed as 

f y | dxdz + f 3 | dxdy - [f y j dxdz] - [f z | dxdy] 

1 y z * y-tAy z+Az 

9T 

~ p dxdydz C D — (2) 

p 9 1 

Applying Taylor’s Series and simplifying, equation (2) becomes. 
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( 3 ) 


Since the principal axes of the material, 5, and ti , are tn general at an angle 
$, with global axes y and z as shown in Pig. 5 , the following transformation 
may be used, 
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we have. 
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( 4 ) 



Substitute (4) into (3), we obtain: 


ORIGINAL pm? r„; 
OF POOR QUAU'ifY 


9 2 T 


(k ? cos 2 {? + k n sin 2 6) 4 - (k ? - k n ) Sin20 


3 2 T 


a^r 


9y 3z 


a 2 t „ ,, 
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Equation (5) is the governing partial differential equation for the heat 
conduction analysis of the blade. 

The boundary conditions around the surface of the blade are: 

(a) prescribed heat flux, q r , in the heat flux input zone, namely, 

q r = 326 BTU/ ft 2 hr . (5; 

(b) convective heat transfer loss, q c , for the entire boundary 
surface of the blade, 

q c = h(T s - T ) (51 

CO 

where T s is the surface temperature of the blade , a variable* 


C* Calculus of Variation 

For a given problem with a governing differential equation and boundary 
conditions, the task of the variational formulation is to find the unknown 
function *F' for the differential equation, which extremizes or makes 
stationary a functional r I* subject to the same boundary conditions of the 


problem. 



where. 


I = f F(x, T(x), 
o 


T' (x)) dx 


T(x), is a weak variation of T(x), namely 


C6) 


T((x) = T(x) + en(x) (7) 

e is an error term of a small magnitude , n(x) is an error function which 
vanishes at its boundaries. 

From the calculus of variation 15], when ’I’ is extremized with respect 
to 'T 1 and the error term e is made to approach zero, the Euler-Lagrange 
equation is formed, namely 


3F d 8F 

— =0 ( 8 ) 

8T dx 3T' 

In other words, the functional ’I’ is extremized or made stationary when 
the Euler-Langrange equation and its boundary conditions are satisfied. 
Therefore, the task of seeking solution, T(x), for the governing partial 
differential equation (5) becomes a numerically simpler problem, that is, 
finding first the ' F' in such a way that its Euler-Lagrange equation is 
identical to the governing differential equation. Next, the functional ’I' 
is formed. Then the finite element method is used to find an approximate 
temperature profile which extremizes the functional and satisfies the 
boundary conditions. The temperature profile so obtained is the solution 
to the heat conduction problem in question. 
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For the qoverninq differential equation, (5), the function ’F 1 may he 
formed as follows: 


1 ST 2 9 T ? 

F = - [{k. cos% + k_ sin 2 fi) ( — ) + {k r sin 2 R + k n cos 2 fl) ( — ) 

2 Q ~ 3y G n 3z 


+ (k^ - k^ ) sin2B 


3T 
( — ) 
3y 


3T 3T 

( — ) + 2p C D — T1 
3z M 3 1 


(9) 


Thus, the functional which contains the function ’F* as described 
in equation (9) and the correspondinq function for boundary conditions 
( 5a ) and { 5b ) , is 


I = f 

v 


FdV + f [q r T + - h {T - T ) 2 ] ds 


(10) 


D. Finite Element Method - variational Approach 

The cross section of the blade at station 1 26 is selected as a 
representative section for temperature profile determination* Trianqular 
elements are used. The blade consists of 3 different kind of materials, 
Douqlas fir, birch plywood and honeycomb paper core. It was divided into 75 
regions, see Fig. 6. Each reqion was further divided into proper numbers of 
columns and rows. See Table 2. Then the trianqular finite element mesh is 
generated automatically [6] • The results of reqions 12 and 54 are shown in 
Fig. 7 as a typical example. 
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To describe the finite element method, a typical element as shown in Pip. R 
is chosen. An approximate temperature profile, T(v,z), can he constructed by 
usinq shape functions N^, and N k to linearly interpolate amonq the nodal 
temperatures T^, T^ and T k , the values of which are to be determined. T(y,z) 
is described ass 


T{y,z) = % Ti + Nj T j + N fc T* 


(11 


where 


Ni = — [aj + b- y + c^ zl 
2A 1 


Ni = — [ai + b ■ y + Ci zl 
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and 


H = - Y k Z j' 
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II 

3~ 
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Substitute T(y,z) into equation (10 ) , and partial differentiate r X 1 
with respect to unknown nodal temperatures. Then, set them equal to 


zero* The result is a set of equations: 



where 
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Note that n is the total number of elements of the blade. [C] is the 
capacitance matrix, [Pj, the heat flux matrix and [K] , the conductance matrix 
of the entire blade. [P e ] and [K e ] contain 'or' in their expressions. The 
purpose is to accommodate the elements with one edge or two edges on the 
boundary of the blade. For example, if edge j-k is on the boundary, the term 
which contains Lj^ is maintained and the other two are deleted. For a non- 
boundary element, j , Lj^ and are assumed a value of aero. Thus, the 
last term of [K e ] is deleted and [P e ] is dropped off for that element. 

Equation (12) describes the heat conduction problem in transient state, 
when it reaches the steady state, the rate of change of temperature respect to 
time vanishes, i.e., 

3{T} _ 

~3t~ “ 


equation ( 12 ) becomes. 


[K] [T] = tP) 


(13) 


equation ( 13 ) describes the heat conduction problem in steady state. 


E. Heat Conduction in Steady State 

The above process is coded into a computer program. To test the program, 
a simple problem on the analysis of one-dimensional heat flow in a rectangular 
slab is tried. The computer solutions agree well to the theoretical solutions 
for both steady state analysis and transient state analysis. 

After the program testing is completed, two steady-state heat conduction 
analyses are made on the University of Toledo’s NAS 6650 computer system. The 
first analysis is performed when the blade section is placed in horizontal 
position and the second, in vertical position. The sun is assumed to shine 
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vertically on the top of the blade, which defines the heat flux input zone as 
shown in Fig, 3 or 4* So, only the element with an edge on the top side has 
heat flux input 'q r ' entering the element thru that edge. The temperatures on 
all nodal points were obtained from the computer solution of equation (13)* 
Nevertheless, only the temperatures of the exterior nodes are plotted in Fig. 

9 and 10. Temperatures of interior nodes are shown representatively on typical 
sections in Fig, 11, 

Observing Fig, 9 and 10, one finds that both of the sunny side surfaces 
reach a maximum temperature of 155°F and the opposite sides are on or slightly 
above 90° F, The identicalness of the results and that the low temperature is 
so close to the ambient temperature offer strong assurances that the analysis 
is correct. Further, Fig. 11 shows a smooth transition from high temperature 
side to low temperature side in all sections. They again indicate that the 
results are quite reasonable. The reliability of the analysis is therefore 
enhanced by engineering judgement. 

F. Heat Conduction in Transient State 

It is important to determine the length of sunshine time required for the 
blade to reach the steady-state temperature. If the length of time required 
is not more than the sunshine hours of a day in the summer, the steady-state 
temperature may be used as the critical temperature distribution for anlaysis. 
Otherwise, the temperature distribution at the end of sunshine hours of a day 
should be used. 

In the time domain, a linear interpolation model is used. The unknown 
temperature field {T}, spanning from one time to another separated by a time 
step At, see Fig. 12, may be expressed as 


{*} = Hi {Ti> + Nj {Tj} 


(14) 
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where and Nj are shape functions, namely. 


»i = 1 


t 

At 


Ni 


At 


hence. 


3Ni 

at" 


-1 

At 


9N j _ 1 

at At 


substitute equation (14) into equation (12) and apply Galerkin's method [7], 
two integral matrix equations are produced. They are, 


At d{T} 

/ N ± UCJ ---- + [K] {T} - {P>) dt = 0 
o 


dt 


(15a) 


At d{T} 

J Nt ( tC] + (K]{T) - {P}) dt = 0 

J dt 


C 15b) 


The time derivative of equation (14) is. 


dN-s 


dN-i 


d{ T} U1, i 

{Ti ) + -d {Tj } 

dt dt dt 3 


- -- {T^ + ~ {T,} 
At At 3 


( 16 ) 


Substitute equations (14) and (16) into equation (15a), we have. 


At t [C] [C] t [K] 

/ (1 )( - {T t } + {T.:} + (1 

0 At At 3 At 1 At 3 

I 

After integration, it yields. 


— ) [KJ {T-;} - {Pi>) dt = 0 
At J 


[C] At [C] At At 

( + — [K] ) {T .) = ( tK] ) {T ± } + — [P ± ] 

26 J i 


2 3 


( 17 ) 
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The constants in parentheses are weighting factors, which may be modified to 
improve the stability of the numerical solution process. Based on Donea's 
suggestion [8], equation (17) is modified to the following form, 


[CJ At 

2 3 


EkI) 



At At 

— [K]) {Ti> + — [P-jJ 
6 2 


08 ) 


Equation (18) is the equation of solution in time domain. 

The initial temperature of the blade is assumed to be uniform and same 
as the ambient temperature, 90°P, namely, 

= {T 0 } = { 90°} 

t=0° 

then, an iterative process is commenced to determine the temperature at the 
next time step t + At, namely {Tj}, by using equation (18). This process is 
repeated until the maximum difference between {Tj} and {T^} is smaller than 
a predetermined small quantity. In this analysis, 500 second is chosen 
as At and good results are obtained as shown in Figs. 13 and 14. 

Observing Figs. 13 and 14 and the computer results, numerical instability 
is experienced in the first hour, but thereafter, stability is restored and 
steady-state condition is reached. The results indicate that the length of 
sunshine time required for the blade to reach the steady-state temperature is 
in the neighborhood of 4 hours. Obviously, this is possible. Thus, steady- 
state temperature distribution of the blade is used as the input of phase II, 
finite element thermal-stress analysis, which will be described in the next 


section 


PHASE IX. THERMAL-STRESS ANALYSIS 


A. Plane-Strain Thermoelastic Problem 

It has been stated earlier that the blade is a slender body, i.e., its 
length is large compared with its maximum cross-sectional dimension. The 
thermal input along the blade axis is close to uniform, namely, the blade 
temperature is independent of the axial coordinate. Under these conditions, 
the use of the concept of plane strain not only provides a good solution [9] , 
but also greatly decreases the computational effort, because it reduces a 
three-dimensional problem to a two-dimensional one. 

If the orthogonal axes system £, n and £ are used, the fundamental 
assumptions for a plane-strain thermoelastic problem are: 

Vl. I = 0, = Ujr(C/Ti)f u n = U n {?fh) (19) 

The corresponding strain components have the form: 

H = Hz = Hr\ = 0 

e CC = ^ ^ ' e nn = e hh^C' 71 E En = s ?ri^' h) (20) 

The stress field is related to the strain field by the elastic constants of 
modulus of elasticity and Poisson's ratios. As stated before, wood is an 
orthotropic material. The modulus of elasticity of wood in the longitudinal 
direction (parallel to grain) is in general, significantly larqer than those 
in radial or tangential directions, in some cases twenty times as large [4], 

To account for the possible difference in radial and tangential directions , a 
10% difference is again maintained in this computation. For plywood, an 
average value is taken. They are listed in Table 3. 


Table 3 » Moflu l\is of Elasticity of Wood 



Douqlas Pir 

Rirch 

Birch Plywood 


1 .95x10 s psi 

2.10x10 s psi 

1.38x10 s psi 

Erp 

0.098X1 0 s pgi 

0.105x10 s psi 

0.77x10 6 psi 

e r 

0.133x10 s psi 

0.164x10 s psi 

0.164x10 s psi 

g tr 

0.014x10 s psi 

0.036x10 S psi 

0.076x10 s psi 


The Poisson's ratios 'v f and the thermal expansion coefficients, a, of wood are 
also siqnificantly different in the major axis directions* They are listed in 
Table 4. 


Table 4* Thermal Expansion Coefficient and Poisson's Ratio of Wood 



Douqlas Pir 

Birch 

Birch Plywood 

a L 

2,1x10”‘ 6 /°F 

2.1x10 S /°F 

1 4.57x1 0~ S 

a T 

34.9X1 0“ 6 /<>f 

39.52x10“ S /°F 

27.05x1 0“ S 

a * 

25.9x1 0** s /op 

30.38x10“ S /°P 

30.38x10~ 6 

V LR 

0.292 

0.426 

0.433 

V LT 

0.449 

0.451 

0*308 

V RT 

0.390 

0.697 

0,476 

V TR 

0.287 

[ 0.447 

0.440 

V RL 

0.020 

0.033 

0.255 

V TL 

0.022 

0.023 

! 0.166 


Reference F 41 * 
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B« Formulation of Potential Energy 

In Phase I, the heat conduction portion of the problem, a qoverninq 
partial differential equation needs to he established first, then, the finite 
element method in variational approach is employed as a numerical procedure to 
obtain its solution# Now, for the problem of findinq the thermal-stresses due 
to a temperature distribution, no control linq differential equation can be 
found. Thus, a conceptually different technique must be used to determine 
directly the displacement field, { f } , without qoinq through the intermediate 
step of establishinq the differential equation# Once f f \ is found, the stress 
field {a } can be determined accordingly# 

First, with an assumed displacement field and under the influence of 
temperature rises at various locations, the total potential enerqy, P.E#, may 
be expressed as 

P.E. = f (- [e l T [El [el - [e] T tE] taTDdv (21) 

v 2 

where V represents the total volume of the blade. In a two-dimensional case, 
it may be replaced by the multiplication of thickness and cross-sectional area, 
i.e. , t * A, 

Let v and w be the displacements of a point in £ and n direction. We may 
express Is] in terms of the displacement field T f > as follows: 

{e }= [0) {f> 


_E c ” 


9 




, tai = 

0 

, ft) = 

V 

E n 


0? 
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-Jen - 


0 

9n 
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9 9 



9n 0? 

- 
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Further, the strain components of a point, fe), may he expressed in terms of 
nodal displacements (D) . For demonstration purpose, a triangular element as 
shown in Fiq. 15 and the aforementioned shape functions, fw] , are again used. 

The strain field, {e} , may he expressed as: 


{e> = [31 m fD>, 
or 

{£} = [B] {D} 

where 

[B] = [91 [Nl 


[NT 


| Ni 0 Nj 0 N k o" 

(O Ni 0 0 N k 


[D] 



v j w i v k 



1 

Ni = — [ai + b£ c + Oi n T 


N a - h ta i “it*')' 1 


"fc “ “ [a k * h k t * °k 1 1 
2A 


Note that A, a^, a-j and a k have been described in (lib) and (11c). 


( 22 ) 


Expanding { D} to cover all the nodal displacements in the entire blade section 
and substituting equation (22) into equation (21), one obtains: 


P.E. 


1 

2 


[D]' r {f[Bl T [E][B]dv> 


[Dl - [n] T r[FO T rEHctT]dv 


( 23 ) 


Equation (23} expresses the total potential enerqv of the blade section in 


finite element formulation 


c 


Finite Element Method — Minimum Potential Energy Approach 


Recall that the total potential energy, p.E*, is formed based on an 
assumed displacement field {f}, for each element. These fields are then 
expressed in terms of unknown nodal displacements by using shape functions# 
Although the displacement fields are admissible, which satisfy the compati- 
bility conditions within the elements and at the nodes, yet the equilibrium 
conditions remain to be established* 

It is known that among all admissible configurations of a conservative 
system, those that satisfy the equations of equilibrium make the potential 
energy stationary with respect to small variations of displacement. If the 
stationary condition is a minimum, the equilibrium state is stable* Thus, 
the equilibrium conditions may be established through the minimization of 
the total potential energy P.E,, namely, 


a p.e* 

= 0 (24) 

3D 

or, (£ / [B] t [E] EBldv) {D} =2 / [B] t [e 3 [<xT]dv (25) 

1 v 1 v 

where v is the volume of an element. 

Equation (25) may be written in the following simplified form: 

[K] [D] = [P] (26) 

where [K] is the total stiffness of the blade and [k] is the stiffness 
of each element, [p e ] is the thermal load produced at nodal points in a 
given temperature field and [E] , Matrix of modulus of elasticity; 


0 RUW& 
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[K1 = E [kl 
1 

(k) = f [B) t [e] [B jdv (27) 


m 

[Pi = Z [p e ] 
1 e 


[p e I = f [B] t [e] [at]dv 
v 

Equation (26) is the finite element formulation of the blade section 
if an isotropic material is used » But the wood composite is an orthotropic 
material, the matrices of modulus of elasticity fEl in equation (25), need 
to be modified accordingly as to be described in the following sections. 


D. The Treatment of Orthotropic Property 

For an orthotropic material such as wood [10] and in the case of plane 
strain, the normal strain-stress relationships apply: 


also 


Thus 



(28) 


( 29 ) 


Note that denotes the strain in the c direction due to strain in 


the n direction 
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Substitution of equation (29) into (20) produces: 

1 1 

e C * - (1 “ ~ <v ?n + v 5n ) o„ 

e n = ” {v nc + v n€ v ee J ff c + ” (1 “ v ne %n } °n 

* h 

Solving the above aquations, one finds: 

a Z = taE ce? + bB ? s n^ f tad “ bc) 

a n = (c Vc + dR n e n> / <ad - be) 

T Sn = G cn Y ?t1 

where, a - 1 - v nC v ?n 

b = v?n + v 5n 

C = V n? + V n g Vpj- 

d = 1 _ v cs v es 

or, the stress-strain relationships may be described in a matrix form as 
(a }= [%] { s } 


where 




1 

=3 

ad - be 


a 

c E, 


b E c 


d P^ 


0 

0 


o 0 (ad - be ) 


From the Maxwell- Betti reciprocal theorem, 

cP h = hF x ' 



so the matrix of modulus of elasticity remains 


( 30 ) 


(31) 


(32) 


(33) 


symmetrical for wood. 
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Equation (33) qives the matrix [T^l for plane strain in an orthotropic 
material, where t; , n and £ are principal directions of orthotropy. 

The elasticity constants of wood, v, E, a and (?, have been described in 
Tables 3 and 4. In view of the blade construction, wood property axis T 
corresponds to axis t, in this section, R to n and L to £ . 

E. Element Axes to Structural Axes 

Because of the blade curvature, the principal axes of material at one 
point in the blade section have in qeneral, a different orientation with 
respect to the structural axes than those at another point. A typical 
orientation of laminated wood has been shown in Fiq. 5, where 5 , n and £ are 
the principal axes of material and local axes of the element, while v, z and x 
are the structural axes and the axes of the blade section. The transformation 
from material axes to structural axes is necessary for all elements in order 
to unify them in the framework of structural axes. This transformation process 
may be described as follows: 

Tjjjg " Transformation from structural axes to material axes 


then, 

{e m l- = [Tjng] {eg) (34) 

where subscript m and s represent material axes and structures axes, and 


cos 2 R 

sin 2 

fl 

sin R 

COS R 

sin 2 8 

cos 2 

B 

-sin R 

cos R 

2 sin 8 cos 

2 sin R 

cos 6 

cos 2 R 

-sin 2 


I 
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In addition, the stress transformation from material axes to structural 
axes is 

! - fw’W t3, > 

l 

t 

f 

f A proper stress-strain transformation may then be obtained by first stating 

j the stress-strain relationships in material axes, namely, 




(36) 


where [E^] has been described in equation (33). Substitution of equation 
(34) and (36) into (35) qives: 

I 

| {*s> = ^ms^ ^ T ms^ ^ G s^ 

I or 


{a g > = [E s ] (e s > (37) 

where 

[E S I = tT ins ] T [E m ] [T ms l 

[E s ] is the matrix of modulus of elasticity for an orthotropic material, of 
which the material axes 5 and n ate oriented at anqle fi with respect to 
principal axes y and a, and material axis 6 coinsides with structural axis x. 
Thisy matrix is used in place of Te] in equation (25) for the determination of 
the stiffness of each element, [k] , and thermal load matrix of each element, 
[P e ]» Then, proper combination of [k] from all elements forms [K], and (P e l 
forms [P] , as described in equations (25) to (27). Finally, equation (26) is 
solved and the displacements of all nodal points [D] are determined. When the 
displacement field is found, the stresses at each element can be computed 
accordingly. Thus, the determination of thermal-stress distribution on the 
wood blade is accomplished. 


F 


Thermal stresses 


n 

f ; The above thermal stress analysis procedure using finite element method, 

is again tested on a simpler problem with known solution. The problem is a 
disk-like structure with a uniformly distributed thermal load applied along a 
diameter. The thermal stresses have been determined by means of mathematical 
; analysis [11], Now, the disk is divided into quadrilateral elements and a 

; computer analysis is performed. The results are in agreement to the mathe- 

matical solution. 

Then, the blade is divided into 388 quadrilateral elements. The 
i nodes are identical to those used in heat conduction analysis, so the nodal 

• i 

temperatures need not be computed again. Corresponding to the two steady-state 
heat conduction analyses made in phase X, two thermal-stress analyses are run 
^ on the University of Toledo computer. The stresses in material axes £ and g 

of all elements are computed. But only the stresses of elements in the surface 
layer are plotted. For example. Fig. 16 shows the stresses in 5 direction of 
the surface layer elements when the blade section is in horizontal position, 

and Fig. 17 shows the stress in the same direction when the blade section is 

in vertical position. 

in order to assure the finite element procedure provides good results, 
the analysis is repeated on a finer-mesh model of 688 elements. For the finer- 
mesh model, the stresses in £ direction of the surface layer elements are 
shown in Fig. 18. A comparison of Fig. 16 to Fig. 18 shows that the stresses 

of the former are converging to the corresponding stresses of the later. This 

is considered to be a good convergence test for finite element analysis. 

Observing Fig. 16 and 17, one finds that on the high temperature sides 
of the blade in horizontal and vertical positions, stresses of same magnitude 
are recorded. This offers an independent check on the correctness of the 








.■r« r* 


computations. The sense of the stresses on the high temperature side is 
tensile, this is largely due to the special characteristics of the blade 
construction. As described before, in D-spar a thin layer of 1/8" birch 
plywood wraps around 3" of lamintated Douglas fir. The thermal expansion 
coefficient in tangential direction of Douglas fir is about 23% higher than 
that of birch plywood. When the temperature rises, the thin layer of birch 
plywood is 'stretched' by the expanding Douglas fir. Thus, tensile stresses 
are produced on the exterior surface of the blade. This result is again in 
agreement with our engineering judgement. Additional rational results that 
give credibility to the analysis are: (1) the smooth transition from high 

stress area to low stress area, (2) in Fig. 16 the stresses on the tail 
portion changing from high to low when the interior material changes from 
Douglas fir to honeycomb paper core; and (3) in Fig. 17 the approximately 
symmetrical stress distribution for a blade of approximately symmetrical 
geometry under an almost symmetrical thermal loading. From the observation 
of the computational results as stated above, one can conclude that, short 
of experimental verification, the analysis offers a set of results with high 


reliability 


27 


Conclusion 

The allowable stresses in principal directions of birch plywood [4] are 
listed in Table 5* A comparison reveals that the maximum blade thermal-stress, 
which occurs in £ direction, is far lower than the allowable stress in the 
corresponding direction, Of course, it is much lower than the ultimate 

strength of the wood. Although the transient thermal-stress analysis is not 
performed due to time limitation of the research project, yet it is believed 
that they should not exceed the allowable stress level for a prolonged period, 
since the sun is a mild heat source, 

Table 5 * Allowable Stresses for Birch Plywood 



Compressive 

Tensile 


5,770 p si 

13,640 p si 

ff T 

3,370 p si 

7,280 p si 

a R 

970 p si 

920 p si 


In conclusion, the cause of cracking of the blade is not due to thermal 
stresses induced by the sun. Investigation of other causes, such as poor 
manufacturing, needs to be made, in addition to the elimination of thermal 
stress as a possible reason of blade cracking, the research brings forth a 
method and presents a practical example in thermal-stress analysis for an 
engineering body of orthotropic materials, which will find many applications 
in the Engineering Community. 
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Appendix 

Computer Programs and Output Examples 

(A) Grid Generation program for Triangular Element 

(B) Bandwidth Reduction program 

(C) Steady State Heat Conduction program - Heat 1 

(D) Transcient State Heat Conduction Program - Heat 2 

(E) Changes on Grid Generation Program, program (A), 

from Triangular Element to Quadrilateral Element 

(F) Changes on Bandwidth Reduction program. Program (B), 

from Triangular Element to Quadrilateral Element 

(G) Transfer Program - preparing Output Data from Programs 

(A) and (C) in the Form of Input Data of SAP IV 

(H) SAP IV Program for Thermal Stress Analysis when the Blade 

Section is in Horizontal position - Quadrilateral Element 

(I) Transformation of Stresses in Structural Axes of SAP IV 

to Stresses in Material Axes 

( j) Output Example 1 : Temperature Distribution when the 

Blade Section is in Horizontal Position 
(K) Output Example 2 s Thermal Stresses of Elements Along 

Boundary When the Blade Section is in Horizontal position 



o n n 


original; pass ts 

OF POOR QUALITY 


Cl 

C| 

C| 

C| 

C| 

C| 

C|‘ 

c 

c 


THIS PROGRAM GENERATES THE GRIDS FOR ANY GIVEN 
REGION 

TRAIKGULAR FINITE ELEMENT 


<> 
<> 
<> 
<> 
<> 

><><><><><> 


DIMENSION NOO( 2500,6) , IRR(500) 

DIMENSION XP(500) ,YP(500) ,XRG(9) , YRG(9) ,N(8) ,NDN(8) ,XC(21 ,21) , 
1JT(74,4) , I COMP (4, 4) ,NN(21 ,21) ,NNRB (74,4,21) , XE (3000 ), YE (3000) , 
1NE(3000), NR(4), LB(3), TITLE ( 10), IIRO(90) ,1100(90), 

2NOU(3000) ,NOUK(3000) ,YC(21 , 21 ) 

REAL N 

DATA NBW/0/ ,NB/0/ , NEL/0/ 

DATA ICOMP/-1.1, 1,-1, 1,-1, -1,1, 1,-1, -1,1, -1,1, 1,-1/ 

DATA IN/5/, 10/6/ 

DATA IP /ll/ 


FINITE ELEMENT GRID GENERATION PROGRAM 


C<><><><><><><><><><><><><><><><><><><><><><><><><><><> <><><><><><><><> 

c 

C TITLE=AN IDENTIFICATION STATEMENT 

C INRG=THE NUMBER OF REGIONS ORIGINALLY DEFINED 

C (CAN SPECIFY A MAX OF 20 REGIONS) 

C INBP=THE NUMBER OF BOUNDARY POINTS OR GLOBAL POINTS ORIGINALLY DEF 

C (MUST HAVE 8 IN EACH REGION) 

C IPCH=0 IF ELEMENT DATA IS NOT TO BE PUNCHED 

C IPCH=1 IF ELEMENT DATA IS TO BE PUNCHED 

C XP(I)=THE X COORDINATE OF EACH OF THE ORIGINAL BOUNDARY POINTS 

C YP(I)=THE Y COORDINATE OF EACH OF THE ORIGINAL BOUNDARY POINTS 

C NRG=THE REGION NUMBER 

C JT(NRG, J)=THE REGION CONNECTIVITY DATA 

C NR0WS=THE NUMBER OF ROWS OF NODES IN A REGION 
C NCOL=THE NUMBER OF COLUMNS OF NODES IN A REGION 

C (CAN SPECIFY A MAX OF 21 ROWS AND 21 COLUMNS IN A GIVEN REGION) 

C NDN(I)=THE GLOBAL NODE NUMBERS ORIGINALLY DEFINED IN A GIVEN REGIO 

C 
C 
C 

100 FORMAT ( 10 A4) 

101 FORMAT (31 3) 

102 FORMAT (8F 6. 2) 

103 FORMAT (513, 2A4) 

1113 FORMAT (IX, 5 110) 

104 F0RMAT(////1X, 10A4//1X, 18HGL0BAL COORDINATES //IX, 

130HNUMBER X COORD Y COORD ) 

105 FORMAT (2X ,I3,7X,F7.4,5X,F7.4) 

106 F0RMAT(//1X, 17HC0NNECTIVITY DATA/ IX, 41HREGI0N SIDE 1 2 

1 3 4 ) 

107 FORMAT(2X,I3 , 14X,4(I2,5X) ) 


A 


o n o 
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108 F0RMAT( 1113) 

109 F0RMAT(///1X,12H*** REGION ,I3,6H ***//10X, 12 ,5H R0WS.10X.I2 

1,8H COLUMNS// 10X, 2 1HB0UNDARY NODE NUMBERS , 1 OX ,815) 

110 FORMATC// IX, 19HREGION NODE NUMBERS/) 

111 FORMAT( IX ,20 15) 

112 FORMAT(//SX, 17HNEL NODE NUMBERS , 9X , 4HX ( 1 ) , 8X,4HY(1) ,8X,4HX(2) ,8X, 

14HY(2) ,8X,4HX(3) ,8X,4HY(3) ) 

113 F0RMAT(1X|5I5,3X,6F12.4) 

114 FQRMAT(4I4 ,SF14 .5) 

115 FORMAT ( / / / IX , 2 1HBANDW IDTH QUANTITY IS,I4,22H CALCULATED IN ELEMENT 
1,14) 

READ (IN ,100) TITLE 

C READ ( IN , 129 )NRA , SEL . RMO , RA , RMI , SHR 

129 F0RMATCI3 , 5F10 .5) 

DO 400 IIM=1 , 1 

NEL=0 

NB=0 

READ (IN, 101) I NRG, INBP,IPCH 
WRITE (10, 116) INRG.INBP.IPCH 

116 FORMAT (IX, 313) 

READ (IN, 102) (XP(I) , 1=1 , INBP) 

WRITE(I0, 102) (XP(I), 1=1, INBP) 

READ(IN , 102) (YP (I), 1=1, INBP) 

WRITE (10, 102) (YP(I), 1=1, INBP) 

DO 1 1=1, I NRG 

READ(IN, 103) NRG, (JT(NRG, J) , J=l,4) ,IRR(I) 

WRITE(IO, 117) NRG, (JTCNRG, J) , J=1 ,4) ,IRR(I) 

117 FORMATC//, IX, 513, 2A4) 

1 CONTINUE 

C WRITE (10, 104) TITLE 

IF(IIM.EQ. 2)G0 TO 401 
WRITE (IP, 100) TITLE 

C WRITE (IP , 129 ) NRA , SEL, RMO , RA , RMI , SHR 

401 WRITE (10, 105) (I ,XP(I) ,YP(I) , 1=1 , INBP) 

C%% WRITE (10, 106) 

DO 2 1=1 , INRG 

C%% WRITE (10, 107) I , ( JT(I , J) , J=1 ,4) 

2 CONTINUE 


V * •''« Jj iW* - L JLA h-L 

C 

C LOOP ON THE REGIONS IN ORDER TO GENERATE THE ELEMENTS 
C 

jV A * t *^ l --V*^* , j*.'**^J**.L*L*.' 1 - J* fc-U -U »* »JL ^ JU « 

** ^ *k ** 4 k 4k *k 4k 4m mm f* 0 % 4\ *k 4k f* r-m »» « rib 4* 4k *% #% 4k 4k 4k 4k 4m- 4% 4\ 4¥ Wk 4k 4k 4k 4m *% 4% 4% 4% Mk 0k 4k 4k 4k 4k 4k 4k 4k 4m. 4m 


DO 22 KK=1,INRG 

READ (IN, 108) NRG, IIRO(KK), IICO(KK) ,(NDN(I) ,1=1,8) 
WRITE (10, 109) NRG, IIRO(KK) ,IICO(KK) , (NDN(I) ,1=1,8) 
NR0WS=IIR0(KK) 

NCOL=IICO(KK) 




GENERATION OF THE ELEMENT NODAL COORDINATES 


A2 


DO 3 1=1,8 


ORK3?!5& 
of pooh o 


I": r, , 


I I=NDN ( I ) 

XRG(I)=XP(II) 

3 YRG(I)=YP(II) 

XRG(9)=XRGC1) 

YRG(9)=YRG(1) 

TR=NR0WS - 1 
DETA=2.0/TR 
TR=NC0L- 1 
DSI=2.0/TR 

DO 4 1=1, NROWS 
TR=I - 1 

ETA=1 . 0 -TR*DETA 
DO 4 J=1 ,NCOL 
TR=J-1 

SI=- 1 , 0+TR*DSI 

N( 1)=- . 25* ( 1 . 0 -SI )* ( 1 . 0 -ETA)* (S I+ETA+1 . 0 ) 
N ( 2 ) = . 5 0 * ( 1 . 0 - S I ** 2 ) * ( 1 . 0 - ETA ) 
N(3)=.25*(1.0+SI)*(1.0-ETA)*(SI -ETA-l.O) 

N (4) = . 50* ( 1 . O+S I ) * ( 1 . 0 -ETA** 2) 

N(5)=. 25*(1 . 0+St )*( 1 . 0+ETA)* (SI+ETA- 1.0) 
N(6)=.50*(l . 0-SI**2)*(l .0+ETA) 

N(7)=. 25*( 1 . 0-SI)*(l . Q+ETA)* (ETA-SI -1 .0) 
N(8)=.S0*(1 . 0-SI )*(1 . 0-ETA**2) 

XC(I,J)=0.0 
YC(I , J)=0 . 0 
DO 4 K=1 , 8 

XC(I,J)=XC(I , J ) +XRG (K)*N (K) 

YC ( I , J)=YC ( I , J)+YRG(K)*N (K) 

4 CONTINUE 

GENERATION OF THE REGION NODE NUMBERS 

KN1=1 

KS1=1 

KN2=NR0WS 

KS2=NCOL 

DO 11 1=1,4 

NRT=JT(NRG,I) 

IF (NRT . EQ . 0 . OR . NRT . GT . NRG ) GO TO 11 
DO 5 J=1 ,4 

5 IF(JT(NRT,J) .EQ.NRG) NRTS=J 
K=NC0L 

IF ( I . EQ . 2 . OR . I . EQ . 4) K=NR0WS 
JL=1 

JK=ICOMP ( I , NRTS ) 

IF(JK.EQ.-l) JL=K 
DO 10 J=1,K 
GO TO (6,7, 8, 9), I 

6 NN (NROWS , J)=NNRB (NRT , NRTS , JL) 

WRITE ( 6 , 7 0 0 0 ) NN ( NROWS , J) 

KN2=NR0WS-1 
GO TO 10 

7 NN(J,NCOL)=NNRB(NRT,NRTS , JL) 

WRITE ( 6 , 7 0 0 1 ) NN ( J , NCO L ) 


A3 


non nan non 


7000 

FORMAT (//' 

STATION 

2 ' 

> 19/) 

7001 

FORMATC//’ 

STATION 

1 ' 

, 19/) 

7002 

FORMATC//' 

STATION 

3 ' 

, 19/) 

7003 

FORMATC//' 

STATION 

4 ' 

, 19/) 


KS2=NC0L-1 
GO TO 10 

8 NN ( 1 , J )=NNRB (NRT ,NRTS , JL) 

C WRITE (6, 7004 )NRT,NRTS,JL 

7004 FORMATC//' NRT=' ,19/' NRTS=\I9/’ 
C VRITE(6,7002)NN(1, J) 

KN1=2 
GO TO 10 

9 NN(J , 1 )=NNRB (NRT , NETS , JL) 

C WR ITE (6 , 7 003 ) NN ( J , 1 ) 

KS1=2 

10 JL=JL+JK 

11 CONTINUE 

IF(KN1 -GT.KN2) GO TO 16 
IFCKS1.GT.KS2) GO TO 16 
DO 12 I=KN1,KN2 
DO 12 J=KS1,KS2 
NB=NB+1 

C WRITE(IO,777)NB 

12 NN(I, J)=NB 

777 FORMATC////' NB = ',16) 
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JL = ' ,19) 


STORAGE OF THE BOUNDARY NODE NUMBERS 

DO 13 I=1,NC0L 

NNRB (NRG , 1, I)=NN (NROWS , I ) 

13 NNRB (NRG ,3,1 )=NN (1,1) 

DO 14 1=1, NROWS 

NNRB (NRG , 2 , I ) =NN ( I , NCOL) 

14 NNRB (NRG , 4 , 1 )=NN(1 , 1) 

OUTPUT OF THE REGION NODE NUMBERS 

WRITE (10, 110) 

DO 15 1=1, NROWS 

WRITE (10,111) (NN ( I , J) , J=1 , NCOL ) 

15 CONTINUE 


DIVISION INTO TRIANGULAR ELEMENTS 

16 CONTINUE 

WRITE (10, 112) 

K=1 

DO 17 1=1, NROWS 
DO 17 J=1,NC0L 
XE(K)=XC(I , J) 

YE(K)=YC(I, J) 

NE(K)=NN(I,J) 

17 K=K+1 
L=NR0WS-1 


DO 21 1=1, L 
DO 21 J=2 , NCOL 

DIAG1=SQRT((XC(I,J)-XC(I+1, J-1))**2+(YC(I, J)-YC(I+1, J-l))**2) 
DIAG2=SQRT( (XC ( 1+1 , J) -XC ( I , J- 1 ) )**2+( YC (1+ 1 , J) - YC f I , J- 1 ) )**2 ) 

NR ( 1 )=NCOL*I+J- 1 
NR (2 )=NCOL* I+J 
NR(3)=NCOL*(I-l)+J 
NR(4)=NC0L*(I-1)+J-1 
DO 21 I J=1 ,2 
NEL=NEL+1 

IF((DIAG1/DIAG2) .GT. 1.02) GO TO 18 
J1=NR(1) 

J2-NR ( I J+l ) 

J3=NR(IJ+2) 

GO TO 19 

18 J1=NR(IJ) 

J2=NR(I J+l) 

J3=NR(4) 

19 LB(1)=IABS(NE(J1)-NE(J2))+1 
LB (2)=IABS (NE (J2)-NE ( J3 ) )+l 
LB (3 )=IABS (NE ( J1 ) -NE (J3 ) )+l 
DO 20 IK=1 ; 3 

IF (LB (IK) .LE.NBW) GO TO 20 
NBW=LB(IK) 

NELBW=NEL 

20 CONTINUE 

IF ( KK . NE . 46 ) GO TO 888 
C WRITE (6, 666) J3 , NE(J3) 

666 FORMAT (///// * J3 = ’,19,' NE(J3) = ’ , 19 ) 

888 CONTINUE 

WRITE (10, 113) KK,NEL,NE(J1) ,NE(J2) ,NE(J3) ,XE(J1) ,YE(J1) ,XE(J2) , 

1 YE(J2) ,XE(J3) ,YE(J3) 

IF(IPCH.EQ.O) GO TO 21 

WRITE(IP, 114) NEL,NE(J1) ,NE(J2) ,NE(J3) ,XE(J1) ,YE(J1) ,XE(J2) ,YE(J2) 
1 ,XE(J3) ,YE(J3) 

NOQ(NEL, 1)=KK 

N00(NEL,2)=NEL 

N00(NEL,3)=NE(J1) 

NOO (NEL, 4)=NE ( J2 ) 

N00(NEL,5)=NE(J3) 

N00(NEL,6)=IRR(KK) 

21 CONTINUE 

22 CONTINUE 

WRITE (10, 115) NBW,NELBW 
DO 95 1B=1,1000 
95 N0U(IB)=0 

C%% WRITE (10, 35) 

35 FORMAT (// / , 10X, * **** BOUNDARY NODE NUMBERS ****') 

40 FORMAT (2213) 

DO 30 1=1 , I NRG 

C%% W r RITE(I0,34)I,IIR0(I) ,IICO(I) 

34 F0RMAT(/3X, ’REGION' ,I2,3X,I2, IX, 'ROWS’ ,3X,I2,1X, ’COLUMNS') 

DO 30 11=1,4 
NRT= JT ( I , I I ) 
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IF(NRT.GT.O)GO TO 30 
K=IIC0(I) 

I1'(II.EQ.2.0R. II.EQ.4)K=IIR0(I) 

C %% WRITE (TO, 33) II , (NNRB(I,II , III) , III“1 ,K) 

33 FORMAT( IX, ' SIDE* , 12 ,/ ,3X, 13 14) 

DO 90 III=1,K 
WI I=NNRB (I, II, III) 

90 NOU(NII)=l 

30 CONTINUE 

IKK=1 

DO 91 IB=1 , 700 
IF(NOUCIB) .EQ.0)G0 TO 91 
NOUK(IKK)=IB 
IKK=IKK+1 

91 CONTINUE 
IKK=IKK- 1 
WRITE (IP, 94) IKK 

94 FORMATC/50X, I4,46X) 

WRITE(IP,92) (NOUK(I) ,1=1, IKK) 

C%% WR ITE ( 1 0 , 9 3 ) (NOUK (I) ,1-1, IKK) 

92 FORMAT (20 15) 

93 FORMAT (15 14) 

DO 333 1=1, NEL 

WRITE(IP,575) ( N00 ( I , J) , J-l , 6 ) 

WRITE(15 ,575) ( NOO(I,J) ,J=1,6) 

C%% WRITE (10, 575) ( N00(I,J) , J=l,6) 

333 CONTINUE 

575 FORMAT ( 5 I 5 , A4 ) 

576 F0RMAT(10X,5I10,10X,A10) 

400 CONTINUE 

STOP 

END 

//* 

/ /* ===== ===== — ■ 1 ■ ===== == ===== 

//* 

//* BANDWIDTH REDUCTION PROGRAM 
//* 

//■ = = = — ========= == 

//* 

c 

COMMON/AAA/ N00(2500,6) ,NTT(80) 

COMMON/ BBB/ X(3, 2000) ,Y(3, 2000) ,NBO(1000) ,TITLE(10) 
COMMON/CCC/ NODES, LMENTS,JT( 12000) ,MEMJT(24000) , JMEM(3000) , 
S JNT ( 3 00 0 ) , IDIFF , MINMAX 
DATA IN/ 11/, 10/6/, IP/ 12/ 

C 

READ (IN, 100)TITLE 
WRITE (IP, 100)TITLE 
100 FORMAT (10A4) 

122 F0RMAT(I3 ,5F10.5) 

J =1 
JJ=0 

150 READ (IN, 17) (JT(3000*(I-1)+J) , 1=1 , 3) , (X(I , J) ,Y(I , J) , 1=1,3) 
JT(9000+J)=0 


: C ! 
C! 
C! 
Ci 
:C! 
CI 


A Co 


OF POOE3 QUAUr 


N0DES=MAX0(JT(J) , JT(3000+J) , JT(6000+J) ,JJ) 

JJ=NODES 

IF(JT(3000+J).EQ,0)GO TO 152 
J=J+1 
GO TO 150 
152 LMENTS=J-1 
NODES=JJ 

12 FORMAT(/ / , 5X , 15HNUMBER OF NODES, 14, 15X, 

1 ' NUMBER OF ELEMENTS ' , 14 , / / ) 

WR;TE(I0,105)TITLE 
105 FORMAT(//,1X,20A4,//) 

WRITE (10 , 12) NODES , LMENTS 
WRITE (10, 13) 

13 FORMAT (// , 3X , 18HNEL JTl JT2 JT3,10X,4HX(1) ,8X, 

$4HY(1) ,8X,4HX(2) ,8X,4HY(2) , 8X,4HX(3) ,8X,4HY(3)) 

17 F0RMAT(4X,3I4,6F14.5) 

DO 10 J=l, LMENTS 

C VRITE(I0,11)J, ( JT(3Q00*( I- 1 )+J) , 1=1,3) , (X(I, J) ,Y(I , J) ,1=1,3) 

11 FORMAT (IX, 415 ,3X,6F12 .4) 

10 CONTINUE 

READ(IN,300)IBO 
300 F0RMAT(50X,I4,46X) 

READ(IN,301) (NBO(I) , 1=1, IBO) 

WRITE (10, 302) 

DO 333 1=1, LMENTS 

READ (IN, 444) (NOQ(I , J) ,J=1 ,6) 

333 CONTINUE 

C WRITE (10, 303 ) (NBO(I) , 1=1 , IBO) 

302 FORMAT (//, * *** BOUNDARY NUMBERS ***' ) 

CALL SETUP 
NTBAN=IDIFF+1 
WRITE (10, 202)NTBAN 

202 FORMAT (// , 2X , 25HTHE ORIGINAL BANDWIDTH IS, 14,//) 

CALL OPTNUM 

IF ( IDIFF . LE . MINMAX) GO TO 115 

MIBAN=MINMAX+1 

lvRITE(IO, 198) 

198 FQRMAT(1H1///,1X,5(3H0LD,3X,3HNEW,7X),/,1X, 

$5 (4HN0DE , 2X , 4HN0DE , 6X) ) 

WRITE (10, 200) (J, JNT(J) ,J=1, NODES) 

WRITE (7, 200) (J, JNT(J) , J=1 , NODES) 

200 FORMAT (5 (15 , IX , 15 , 5X) ) 

WRITE (10, 13) 

DO 180 J=l, LMENTS 
DO 555 11=1,3 
I I 1=11+2 

555 NOO(J, III)=JNT(JT(3000*(II-1)+J) ) 

C 

180 CONTINUE 

C WRITE (10, 11) J, (JNT(JT(30QQ*(I-1)+J)) , 1=1,3) , 

C $(X(I, J) ,Y(I, J) ,1=1,3) 

WRITE (10, 201)MIBAN 

201 F0RMAT(//,2X, 'THE NEW BANDWIDTH IS', 13,/) 

WRITE ( IP , 3 10)NODES , LMENTS , MI BAN 


A7 


DO 181 J=l, LMENTS 0P Qc!i-U-*j “ 

181 WRITE(IP, 18) J, (JNT(JT(3000*(I-1)+J)) , 1=1,3) , (X(I, J) , 
$Y(I,J),I=1,3) 

18 FORMAT(4I5,6F10.4) 

WRITE (10 ,302) 

WRITE (IP , 300 ) IBO 

WRITE (IP, 301) (JNT(NB0(I) ) , 1=1 , IBO) 

C WRITE (10 ,303) (JNT(NBO(I ) ) , 1=1 , IBO) 

303 F0RMAT(15I4) 

GO TO 117 
115 CONTINUE 

WRITE (10, 101) 

117 CONTINUE 

101 F0RMAT(///,2X,'THE ORIGINAL BANDWIDTH IS MIN'INMUM* ) 

WHITE ( IP , 3 10) NODES , LMENTS , NTBAN 
310 FORMAT (3 15) 

DO 190 J=l, LMENTS 

190 WRITE (IP , 18) J, (JT(3000*(I-1 )+ J) ,1=1,3) , 

1(X(I , J) ,Y(I , J) ,1=1,3) 

WRITE(IP,300)IB0 

WRITE(IP,301)(NB0(I),I=1,IB0) 

301 FORMAT (2015) 

C C! 

444 F0RMAT(5I5 ,A4) 

LLL = NOO (1,6) 

C 

L0L=0 

MMM=0 

NNN=1 

DO 666 I 1=1, LMENTS 

IF ( LLL. EQ. NOO (II ,6) )G0 TO 1888 

NNN=NNN-1 

LOL=LOL+NNN 

MMM=MMM+1 

NTT (MMM)=LOL 

C WHITE ( 10 , 8 8 8 ) NNN , NTT ( MMM ) 

LLL=NOO(II ,6) 

NNN=1 

1888 CONTINUE 

C WRITE(IO,777)NNN, (NOO(II , J) , J=1 ,6) , (X(1 ,11) ,Y(I , II) ,1=1,3) 

IF (II. NE . LMENTS ) GO TO 1889 
LOL=LOL+NNN 
MMM=MMM+1 
NTT(MMM)=LOL 

C WHITE (10, 888) NNN,NTT(MMM) 

1889 NNN=NNN+1 
666 CONTINUE 

C 

777 FORMAT(4X,6I6,4X,A8,6F12.5) 

888 FORMAT(//’ NUMBER OF ELEMENTS IN THIS REGION =’,2110/) 

400 CONTINUE 
C 

C - 

C - - 


AS 


omcmiL 


OF POOH: QoALi 


STOP 

END 

SUBROUTINE OPTNUM 
C 

DIMENSION NEWJTC3000) , JOINT(3000) 

COMMON/AAA/ N00(1500,6) ,NTT(80) 

COMMON/BBB/ X(3,2000) ,Y(3,2000) ,NB0(1000) ,TITLE(10) 

COMMON/ CCC/ NODES , LMENTS , JT ( 12000) ,MEMJT(24000) , JMEM{3000) , 
$JNT(3000) , ID IFF .MIN'MAX 
MINMAX=IDIFF 
DO 60 IK=1, NODES 
DO 20 J =1 .NODES 
JOINT (J)=0 
20 NEWJT(J)=0 

MAX=0 
1=1 

NElvJT ( 1 )=IK 
JOINTf IK)=1 
K=1 

30 K4=JMEM (NEWJT ( I ) ) 

IF(K4.EQ. 0)G0 TO 45 
JSUB= (NEW JT ( I ) - 1 ) * 8 
DO 40 JJ=1 ,K4 
K5=MEMJT(JSUB+JJ) 

IF(JOINT(K5) . GT. 0)00 TO 40 
K=K+1 

NEWJT(K)=K5 
JOINT CK5)=K 
NDIFF=IABS(I-K) 

IF(NDIFF.GE.MINMAX) GO TO 60 
IF (NDIFF . GT . MAX) MAX=NDIFF 
40 CONTINUE 

IF (K.EQ. NODES) GO TO 50 
45 1=1+1 

GO TO 30 

50 MINMAX=MAX 

DO 55 J=l, NODES 
55 JNT(J)=JOINT(J) 

60 CONTINUE 

RETURN 
END 

SUBROUTINE SETUP 

COMMON/ AAA/ N00(1500,6) .NTT (80) 

COMMON/BBB/ X(3.200G) ,Y(3 ,2000) .NBO(IOOO) .TITLE (10) 
COMMON/CCC/ NODES, LMENTS . JT(12000) ,MEMJT(24000) , JMEM(3000) , 
?JNT(3000) , IDIFF.MINMAX 
IDIFF=0 

DO 15 J=l, NODES 
15 JMEM(J)=0 

DO 60 J=l, LMENTS 
DO 50 1=1,4 
JNTI=JT(3000*(I-1)+J) 

IF(JNTI.EQ.O) GO TO 60 
JSUB=( JNTI - 1 ) *8 
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fto 


20 

30 


40 

50 

60 


C 

C 

C 

C 

C 

C&1B 

C 


ORIGINAL PPf;:£ IT 

DO 40 11=1,4 0F POOR QUAiiiY 

IF(II.EQ.I) GO TO 40 

JJT=JT(3000*(II-1)+J) 

IF(JJT.EQ.O) GO TO 50 
MEM 1=JMEH ( JNTI ) 

IF(MEMl.EQ.O) GO TO 30 
DO 20 1 1 1=1 , MEM 1 

IF(MEMJT(JSUB-fIII).EQ.JJT) GO TO 40 
CONTINUE 

JMEM C JNTI ) =JMEM ( JNTI )+l 
MEMJT ( J SUB+ JMEM ( JNTI ) ) = J JT 

IF C IABS C JNTI - JJT) . GT . IDIFF) I DIFF= IABS ( JNTI - JJT) 

CONTINUE 

CONTINUE 

CONTINUE 


RETURN 


END 




' f l IjJ* *'■ ■ 


HEAT . . . 1 


JL +.** 


<> 

<> 


c <> 

REAL KXX,KXXX 

DIMENSION TITLE (10) ,NT1(11) ,NT2(80) 

DIMENSION NS(3),T(3) 

DIMENSION IS IDE (2) ,IC0N(2) ,PHI (3) 

DIMENSION CFV ( 6 ) , NODNUM ( 6 ) 


INTEGER NP,IER,N 

REAL*8 Z(482) ,R(482) ,AREAE(752) ,X(3) ,Y(3) , DDT, TIME, TDT 
2 ,X1,X2,X3 , Y1 , Y2,Y3 ,C0F,THSTF(482 , 35) ,THCAP(482 , 35) 
COMMON/CONVEN/H ( 752) , ISIDEH (752,2; 

C 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 ,THCAP,NPT(752,3) 

REAL*8 ZZ(482,35),WK(8694),F,EF 
COMMON / SOLVEQ / F(482, 1) ,EF(752) 

REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /CONT/ KZETA '752) ,KETA(752) ,RH0CP(752) ,BETA(752) 

C 

DATA IN/ 12/, 10/6/ 

C DATA NT1/343, 347, 351, 355, 359, 363, 367, 371, 375, 379, 383/ 

C 

C DEF'NITION OF THE CONTROL PARAMETERS 

C AREAE ( NE ) : ARRAY CONTAINS AREA OF EACH ELEMENTS 

C R ( NE*3): ARRAY CONTAINS THE X-COORDINATES OF EACH ELEMENT 

C Z ( NE*3) : ARRAY CONTAINS THE Y-COORDINATES OF EACH ELEMENT 

C T ( 3 ) : ARRAY CONTAINS THE TEMP. OF THE ELEMENT'S THREE 


Ai o 
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C NODES . 

C TOLD ( NE*3) : ARRAY COPIES THE T(3) ARRAY FOR EACH ELEMENT 

C NE = NUMBER OF ELEMENTS 

C NP = NUMBER OF GLOBAL TEMPERATURES 

C NBW = BAND WIDTH 

C H = CONVECTION COEFFICIENT 

C TINT = AMBIENT AIR TEMP. 

C 

C 

c 

c 

C INPUT OF THE TITLE CARD AND THE CONTROL PARAMETERS 
C 

TINF =10.00 
REWIND 12 

C 

READ (IN'. 3) TITLE 

C WRITE (10,4 ) TITLE 

READ ( IN , 2 ) NP , NE , NBW 

NBW=NBW+17 

WRITE (10, 11021 NP , NE , NBW 

C 

1 F0RMAT(4I5 ,6F10.4) 

2 FORMAT (3 I 5 ) 

3 FORMAT (10A4) 

4 FORMAT (1H1,///1X,1 0A4 ) 

1102 F0RMAT(1X, 316) 

C 

C 

C 

C OUTPUT OF TITLE AND DATA HEADINGS 

C - 

44 FORMAT ( 1 HI , // / / IX , 20 A4/ / IX , 5HKXX =,F7 . 1, 10X.5HKYY = .F7.1/1X, 
121HCONVECTION COEFF(l) =,F7 . 1 , 3X,21HC0NVECTI0N C0EFF(2) =,F7.1/ 
21X, 15HFLUID TEMP(l) =,F7 . 1 ,3X, 15HFLUID TEMP (2) =,F7.1//,1X, 
317HNEL NODE NUMBER, 6X,4HX(1) ,6X,4HY(1) ,6X,4HX(2) ,6X,4HY(2) , 
46X,4HX(3),6X,4HY(3) ) 

C ASSEMBLYING OF THE GLOBAL STIFFNESS MATRIX AND GLOBAL FORCE MATRIX 

C********* 

c 

C INPUT AND ECHO PRINT OF ELEMENT DATA 
C 

C 

C WRITE (6, 604) 

604 FORMAT (1H1, /////) 

C WRITE (6, 304) 

304 FORMAT ( * T(l) T(2) T(3) 

1 ISIDE(l) ISIDE(2) ’/) 

1061 FORMAT (6F 10. 3, 21 10) 

11 FORMATC ' ,3F10.3,2I10) 

WRITE (6 ,1107) 

WRITE (6 , 1108) 


C 

C 


ORIGINAL P/V T : r w 
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c 

C NE : TOTAL NUMBER OF ELEMENTS 

C 

C. .IN: 12 MAIN LOOP 

T(l)=10. 

T(2)=10. 

T(3)=10. 

DO ISO 1=1, NE 
H(I)=0 .001D0 
ISIDEHd , 1)=0 
ISIDEH(I,2)=0 
180 CONTINUE 

DO 178 JJ=1 ,90 
N=2“ JJ 

ISIDEHfN, 1)=2 
ISIDEH(N,2)=0 
H(N)=5 .DO 
178 CONTINUE 
C 

C BIRCH LOOP 

C 

DO 1001 J=1,2B0 
KZETA(J)=. 0924 
KETA(J)=. 0836 

1001 CONTINUE 
C 

C DOUGLAS FIR LOOP 

C 

DO 1002 J=281, 302 
KZETA(J)=. 0754 
KETA(J)=. 06820 

1002 CONTINUE 

DO 1010 J=35 1,638 
KZETACJ)=.07S4 
KETA(J)=. 0682 
1010 CONTINUE 

DO 1015 J=687 , 752 
KZETA(J)=.0754 
KETA(J)=. 0682 
1015 CONTINUE 
C 

C PAPER LOOP 

C 

DO 1020 J=303 ,350 
KZETA(J)=. 011161 
KETA(J)=. 011161 
1020 CONTINUE 

DO 1111 J =639,686 
KZETA(J)=. 011161 
KETA(J)=. 011161 
1111 CONTINUE 
C 

C LOOP ON THE REGIONS 

C 


Al2 


n o 


O ;• 
OP POOR 

DO 7 KK=1 ,NE 

IF(KK.GE. 1. AND. KK.LE. 22) BETA (KK)=- .2618 
IF (KK . GE . 23 . AND . KK . LE . 34) BETA (KK)=- . 15 707 
IF (KK . GE . 35 . AND . KK . LE . 58) BETA(KK)=0 . 0 
IF (KK . GE . 59 . AND . KK . LE . 66) BETA (KK)=. 38394 
IF (KK . GE . 6 7 . AND . KK . LE . 74) BETA(KK)=. 17452 
IF(KK. GE. 75 . AND.KK. LE ,82) BETA (KK) =1 . 22164 
IF(KK . GE . 83 . AND . KK . LE . 90) BETA(KK)=. 69808 
IF (KK . GE . 9 1 . AND . KK . LE . 98) BETA(KK)=2 . 4433 
IF(KK.GE.99.AND.KK.LE. 106) BETA(KK)=1 . 91972 
IF(KK.GE . 107 . AND. KK.LE .114) BETA (KK) =2. 9668 
IF (KK . GE . 1 15 . AND . KK . LE . 122) 

IFCKK.GE. 123. AND. KK.LE. 146) 

IF(KK.GE. 147. AND. KK.LE. 158) 

IF (KK.GE. 159. AND. KK.LE. 202) 

IF(KK.GE. 203. AND. KK.LE. 214) 

IF(KK.GE. 215. AND. KK.LE. 230) 

IF (KK . GE . 23 1 . AND . KK . LE . 246) 

IF (KK.GE. 247. AND. KK.LE. 258) 

IF (KK . GE . 25 9 . AND . KK . LE . 326) 

IF(KK.GE. 327. AND. KK.LE. 350) BETA (KK) =-. 15707 
IFCKK.GE. 351. AND. KK.LE. 430) BETA(KK)=0.Q 
IF(KK.GE. 431. AND. KK.LE. 438) 

IF (KK . GE . 439 . AND . KK . LE . 446 ) 

IFCKK.GE. 447. AND. KK.LE. 454) BETA(KK)=. 38394 
IF(KK . GE . 455 . AND . KK . LE . 462) BETA (KK)= . 17452 
IFCKK.GE. 463. AND. KK.LE. 470) BETA(KK)=1 . 2216 
IFCKK.GE. 471. AND. KK.LE. 478) BETA (KK)=. 69808 
IFCKK.GE. 479. AND. KK.LE. 486) BETA(KK)=1.2216 
IFCKK.GE. 487. AND. KK.LE. 494) BETA(KK)=. 69808 
IF (KK.GE. 495. AND. KK.LE. 502) 

IFCKK.GE. 503. AND. KK.LE. 510) 

IFCKK.GE. 511. AND. KK.LE. 518) 

IF (KK.GE. 519. AND. KK.LE. 526) 

IFCKK.GE. 527. AND. KK.LE. 534) 

IFCKK.GE. 535. AND. KK.LE. 542) 

IF (KK.GE. 543. AND. KK.LE. 550) 

IF (KK.GE. 551. AND. KK.LE. 558) 

IF (KK.GE. 559. AND. KK.LE. 638) 

IFCKK.GE. 639. AND. KK.LE. 662) 

IF (KK . GE . 663 . AND . KK . LE . 708) BETA(KK)=3 .40314 
IF (KK.GE. 709. AND. KK.LE. 752) BETA(KK)=0.0 


BETA (KK) =2. 7574 
BETA(KK)=3. 14136 
BETA (KK)=3. 29843 
BETA(KK)=3 .40314 
BETA (KK)=3. 29843 
BETA(KK)=3. 14136 
BETA (KK) =0.0 
BETA(KK)=- . 15707 
BETA (KK}=- .2618 


BETA(KK)=. 38394 
BETA (KK)=. 17452 


BETA(KK)=2 . 4433 
BETA(KK)=1.91972 
BETA (KK)=2, 4433 
BETA(KK)=1. 91972 
BETA (KK)=2. 9668 
BETA (KK)=2. 7574 
BETA (KK)=2. 9668 
BETA(KK)=2 . 7574 
BETA(KK)=3. 14136 
BETA (KK)=3. 29843 


121 

122 


READ(IN, 1) NEL,NS,X(1) ,Y(1) >X(2) ,Y(2),X(3),Y(3) 

WRITE (6, 121) NEL,NS,X(1) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) ,BETA(KK) 
WRITE (7, 122) NEL,N5,X(1) ,Y(1) ,X(2) ,Y(2) ,X(3) , Y(3) , BETA(KK) 
FORMAT (415, 7F10. 4) 

FORMAT (414, 7F8. 4) 


DIMENSION IN FEET... 
X1=X(1) 

X2=X(2) 

X3=X(3) 

Y1=Y ( 1 ) 
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Y2=Y(2) 

Y3=Y(3) 



X 

C 

1 

)= 

X(l) 

/12.OD0 

X 

( 

2 

)= 

X(2) 

/12.0D0 

X 

( 

3 

)= 

X(3) 

/12.0D0 

Y 

( 

1 

)= 

Y(l) 

/12.0D0 

Y 

( 

2 

)= 

Y(2) 

/12.0D0 

Y 

( 

3 

)= 

Y(3) 

/12.0D0 
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T1=T(1) 

T2=T(2) 

T3=T(3) 

1107 FORMAT (1H1 ,///20X 

* , ' T(l) T(2) T(3) ' 

* ISIDE(l) ISIDE(2) ’) 

1108 FORMATC* * , 20X 

* 1 - 

> 

j. i i ^ 

C WRITE(I0,n)T,ISIDEH(KK,l),ISIDEHCKK,2) 

c 

c 

c 

C T0LD(NS(1))-T(1) 

C TOLD (NS (2 ) )=T( 2 ) 

C TOLD (NS(3))=T(3) 

C H : CONVECTION COEFFICIENT 

C WRITE(6,5152)KK,NS(1) ,NS(2) ,NS(3) ,TOLD(NS(l) ) ,T0LD(NS(2) ) , 

C 1 T0LD(NS(3)) , H(KK) , (ISIDE (JJ) , Jj=l ,2) , 

C 1 Xl,Yl,X2,Y2,X3,Y3,NOT(KK,l),NOTCKK,6) 

5152 FORMATC* * ,415 , 3F10. 3 ,F5 . 0 , 212 , 6F8 A , 15 ,A8) 


C 


T ( 1 )=T1 
T(2)=T2 
T(3)=T3 


ECHO THE INPUT 

WRITE (10, 23) NEL,NS ,X( 1) , Y(l) ,X(2),Y(2) ,X(3) ,Y(3) 

* , HH,RH0, CP, ISIDE ,T 

23 FORMAT(1HO, 15 , IX, 315 , 2X,6(1X, F8 . 4) , 2X, 2F5 .0,F6 . 3, 2X, 213, 

* 2X,3F6 . 0) 

C CALCULATION OF THE CONDUCTION MATRIX 

C 

C NPT (NE , 3 ) : ARRAY STORES THE NODE NUMBERS FOR EACH ELEMENT 
NPT(KK, 1) = NS ( 1 ) 

NPT(KK,2) = KS(2) 

NPT(KK,3) = NS ( 3 ) 

RCNS(l)) = XC1) 

R(NS(2)) = X(2) 

R(NS(3)) = X(3) 

ZCNS(l)) =Y(1) 

Z(NS(2)) =Y(2) 

Z(NS(3)) =Y(3) 

AREAE(KK)=CXC2)*Y(3)-X(2)*YC1)+X(3)*Y(1)-XC3)*Y(2) 


ORIGINAL PASS 
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1 + X(l)*Y(2)~X(l)*Y(3))/2.0 


C i! 

C 

C WRITE (6 ,266)KK,X(1) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) .AREAE(KK) 

266 F0RMAT(1H0,I5 ,6(1X,F8 .4 ),2X,D14.7) 

C WRITE (6, 3 19 )KK, (T(I), 1=1,3) 

319 FORMAT! ' ELEMENT NO. \I8,3F12.4) 

7 CONTINUE 

C END OF THE DO LOOP Cl 

C WRITE (6 , 264) (I ,TOLD(I) ,1=1 ,NP) 

264 FORMAT (1HO ,I5,E14.5,2X,I5,E14.5,2X,I5,E14.5,2X,I5,E14.5 


* , 2X, 15 ,E14 .5 ,2X , 15 ,E14 . 5) 

C 

C ADD EXTERNALLY APPL. CONC. NODAL FORCES VECTOR TO EF VECTOR 
ID1=0 
INK=0 

202 READ (IN, 203) NODNUM.CFV 

203 FORMAT(6I3 , _X,6F10 .5) 

ID=0 

DO 204 L=1 , 6 

IF( NODNUM(L).LE.O) GOTO 205 
ID=ID+1 

INODNU=NODNUM(L) 

204 EF ( INODNU ) =CF V ( L ) +E F ( INODNU ) 

GOTO 206 

205 INK=1 
IF(ID.EQ.O) GOTO 216 

206 IF(ID1.EQ.1)G0T0 222 
WRITE (10, 208) 

208 F0RMATC1H1, 1 NODE NUMBER ',10X, ' GENERALIZED NODAL FORCES 

222 WRITE (10, 207) (NODNUM(L) ,CFV(L) ,L=1,ID) 

207 FORMAT( 1H0 , 6 (13 , E 14 . 5 , 2X) ) 

IF ( INK . EQ . 1 ) GOTO 216 
ID1=1 

GOTO 202 
216 CONTINUE 

C 

JJJKNT=0 

C 

CALL FORMTK (NE.NP.NBW) 

CALL FORMTF (NE ,NP,NBW) 

CALL ASSMBL (NP,NBW,ZZ) 

C WRITE(6, 9999) ((ZZ(I, J),J=1,NBW), 1=1, NP) 

9999 F0RMAT(2X, 7G11 .5) 

C WRITE (6,9988) (F(I , 1) , 1=1 ,NP) 

9988 FORMAT (2X,G1 1.5) 

CALL LEQT1B (ZZ , NP , 17 , 17 , NP , F , 1 , NP , 0 , WK , IER ) 

JJJKNT=JJJKNT+1 

IF (JJJKNT . GE . 1 ) GOTO 979 

GOTO 980 

979 CONTINUE 

980 CONTINUE 
WRITE (10 ,666) 


C! 

-C\ 


’) 
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666 FORMAT ( ' NODE NO. ' ,7X, * TEMPERATURE IN DEGREE F* ,/) 

WRITE (6, 787) (I,F(I,1),I=1,NP) 

787 FORMAT ( 2X,I5 , 10X.F15 . 8) 

WRITE (7 , 777) (Ff I , 1) , 1=1 ,NP) 

777 FORMAT(9F8.4) 

C 

C 

STOP 

END 

SUBROUTINE FORMTK (NE.NP.NBW) 

C zWrie “V iV vr i r* & V? VrtV v -V -,V sV * ft ft * ft ft ftft ft ft V,* ft ftftft ft -.V ft ft ft ft ft V: ft -.V -.V ft ft ft -V ft i 

C FORM BANDED THERMAL CONDUCTANCE MATRIX 

C 

IMPLICIT REAL*4(A-H,0-Z) 

REAL*8 KZETA , KETA , RHOCP , BETA 

REAL*8 B I , B J , BK , C I , C J , CK , COF , CON 1 , CON2 , CON3 

REAL*8 RB AR , CON , LG , RCTL , ELKT (3,3), THSTF(4S2 , 35 ) 

REAL* 8 AOO,A11,A22,BOO,B11 ,B22,C00,THCAP(482,35) 

REAL*8 Z (482) ,R(482) ,AREAE(752) ,X(3) ,Y(3) , DDT, TIME ,ATIM,P 
COMMON/ CONVEN/H ( 75 2 ) , I S IDEH ( 75 2 , 2 ) 

C 


COMMON / SOLID / Z ,R,AREAE ,THSTF 
1 ,THCAP,NPT(752 9 3) 

REAL*8 ZZ(482,35) ,F,B,FF,EF 
COMMON / SOLVEQ / F(482 , 1) ,EF(752) 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHOCP(752) , BETAC752) 

C 

C 

C 

C INITIALIZE GLOBAL THERMAL STIFFNESS MATRIX 

. IN=5 
10=6 

DO 120 J=1 , NBW 
DO 110 1=1 ,NP 
THSTF(I.J) = 0 .DO 
110 CONTINUE 
120 CONTINUE 
C 
C 

130 NBAND = 0 

DO 899 J=1,NE 

C WRITE (6 , 788 )KZETA (J) , KETA (J) , BETA (J) 

788 FORMAT ( 3X , 3F20 . 9 ) 

899 CONTINUE 
C 

C WRITE (6, 7776) 

7776 FORMAT ( ’ CONDUCTION MATRIX ELEMENTS ARE ' ,/) 

DO 230 N=1 ,NE 

C<><><><><><> <><><><><><><><><><><> 


II = NPT(N, 1) 

JJ = NPTCN',2) 

KK = NPT(N,3) 

BI = Z(JJ)-Z(KK) 


A 
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555 


Oft’G'M. V' 


BJ = ZCKK)-Z(II) OF POOR Q ,J Ai- - 

BK = Z(II)-Z(JJ) 

Cl = R(KK)-R(JJ) 

CJ = R(II)-R(KK) 

CK = R(JJ)-RCII) 

X(l) =R C 1 1 ) 

X c 2 ) =R( JJ) 

X(3) =R(KK) 

Y(ll “Z(II) 

Y(2) =Z( JJ) 

Y(3) =Z(KK) 

WRITE(6,555)X(1),X(2),X(3),Y(1),Y(2),Y(3) 

FORMAT (2X.6F12.6) 

ALPHA=BETA(N) 

CON 1= ( KZETA (N ) * ( COS C BETA ( N ) ) ) •“• 2+KETA ( N ) * ( S I N ( BETA (N ) ) ) **2 ) 
1/ (4 ,*AREAE(N) ) 


CON2=CKZETA(N)-KETA(N) ) *S IN ( 2 . * B ETA ( N ) ) / (4.*AREAE(N) ) 

CON3= (KZETA ( N ) * ( S I N ( BETA (N ) ) )** 2+KETA ( N ) * ( COS ( BETA c N ) ) ) **2 ) 
1/ (4.*AREA£ (N) ) 

CON 1= (KZETA (N ) * ( COS ( ALPHA ) ) **2+KETA ( N ) * ( S I N (ALPHA ) ) **2 ) 

1/ (4.*AREAE(N) ) 

CON2= (KZETA (N) -KETA (N) )*SIN (2 . *ALPHA) / (4 . *AREAE (N) ) 

CON3= (KZETA (N)* (SIN (ALPHA) )**2+KETA (N)* (COS (ALPHA) )**2 ) 

1/ (4 ,*AREAE (N) ) 

C 

C WRITE (6 , 7 9 9 9 ) CONI , CON2 , CON3 , AREAE (N) 

7999 F0RMAT(2X,3F10. 3, 2X,F10. 7) 

ELKT(1, 1) = ( B 1**2 ) *CON 1+ ( C 1**2 ) *CON3+ ( B I*CI ) *C0N2 
ELKT(1,2) = (BI*BJ)*CONl+(CI*CJ)*CON3+(BI*CJ)*CON2 
ELKT(1,3) = (BI*BK)*C0N1+(CI*CK)*C0N3+(BI*CK)*C0N2 
ELKT (2,1) = (BI*BJ)*CONl+(CJ*CI)*CON3+(BJ*CI)*CON2 
ELKT (2,2) = ( B J**2 ) *CON 1+ ( C J**2 ) *CON3+ ( B J*C J ) *CON2 
ELKT(2,3) = ( B J*BK ) * CON 1+ (C J*CK ) *CON3+ (B J*CK) *CON2 
ELKT (3, 1) = (BK*BI)*CONl-f(CK*CI)*CON3+(BK*CI)*CON2 
ELKT(3 ,2) « ( BK*B J ) *CON 1+ ( CK*C J ) *CON3+ ( BK*C J ) *C ON2 
ELKT (3 ,3) = (BK**2)*CONl+ ( CK— 2 ) --CON3+ (BK*CK)*C0N2 
WRITE(6,*)N,H(N) , ( ISIDEH(N, Jl) , Jl=l,2) 

DO 10 IQ=1,2 

C WRITE (6 , 7805)N, ( ISIDEH(N, Jl) , Jl=l ,2) 

7805 FORMAT ( ' SIDE EXPOSED ’,3110) 

IF ( ISIDEH(NMQ) . LE . 0 ) GO TO 240 
JQ = ISIDEH (N,IQ) 

C WRITE (10 ,12) JQ,N 

12 FORMAT ( 1H0 , ' INSIDE 13 ,3X, ' OF ELEMENT' , 15) 

KQ = JQ +1 

IF ( JQ.EQ . 3) KQ=1 

LG = DSQRT( (X(KQ) -X(JQ) )**2+(Y(KQ) -Y ( JQ) )**2) 

HL =H(N)*LG/6. 

WRITE (6,*) N,H(N) ,LG,HL 
IF(JQ.EQ.l) GO TO 20 
IF (JQ.EQ. 2) GO TO 30 
IF(JQ.EQ.3) GO TO 40 
GO TO 10 

20 ELKT(1 , 1)= ELKT(1,1)+HL * 2. 
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ELKT(2,2)= ELKT(2,2)+HL * 2. 

ELKT(1 ,2)= ELKTCl, 2)+HL 
ELKT(2, 1)=ELKT(2, 1)+HL 
GO TO 10 

30 ELKT(2,2)= ELKT(2,2)+HL*2. 

ELKT(3,3)= ELKT (3,3) +HL* 2 . 

ELKT(2,3)= ELKT(2 ,3)+HL 
ELKT(3,2)=ELKT(3,2)+HL 
GO TO 10 

40 ELKTCl, 1)= ELKTCl, 1)+HL * 2. 

ELKT(3,3)= ELKT(3,3)+HL * 2. 

ELKTCl ,3)= ELKTCl, 3)+HL 
ELKTC3, 1)=ELKT(3, 1)+HL 
10 CONTINUE 

WRITE (10, 2000 )N, II , JJ,KK,X(1) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) ,C0NDTY(S) , 
1AREAE(N),RBAR,C0N,ELKT(1,1),ELKT(1,2),ELKT(1,3),ELKT(2,2), 
2ELKT(2,3),ELKT(3,3) 

2000 FORMATC 1H0 ,415 ,8G11.4,/21X,8G11 .4) 

C J 

STORE IN GLOBAL MATRIX C! 

C! 

240 CONTINUE 

WRITE ( I0,266)N J X(1),YC1),X(2),Y(2),XC3),YC3),AREAE(N),NBAND 
266 FORMATC 1H0 , IS , 6 C IX ,F8 .4) ,2X,D14. 7 ,2X,I5) 

WRITE (6, 27 7 )N, CC ELKTCl, J) , J=l,3) ,1=1,3) 

277 FORMATC' ’,15/ 3(3F25.6/)) 

DO 220 LL=1 , 3 
I=NPT(N,LL) 

DO 210 MM=1. , 3 
J=NPT (N , MM) - 1+18 

THSTF ( I , J) =THSTF ( I , J) +E LKT ( LL , MM ) 

WRITE (6 ,7543)1 ,J, II ,JH,NPTCN,J) ,NPT(N, I) 

7543 FORMATC 1 1 , 618 , ' ****** ' ) 

IF (J.GT.NBAND) NBAND=J 
210 CONTINUE 
220 CONTINUE 
230 CONTINUE 

WRITE (6 , 5050) ( (THSTF (I ,J) , J=1,NBW) , 1=15 , 20) 

5050 FORMATC 7G11. 5) 

WRITE (6,200) NBAND 
200 FORMATC 1H0,' NBAND= ’,15) 

END OF PROGRAM 


RETURN 

END 

SUBROUTINE FORMTF (NE ,NP , NBW) 
FORM THERMAL FORCE VECTOR 
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REAL*8 ELF (3 ) , LG , RCTL , COF , THCAF C482 , 35 ) , THSTF (482 , 35 ) 
REAL*8 Z(482) ,R(482) ,AREAE(752) ,X(3) ,Y(3) ,DDT 
COMMON/CONVEN/H (752), ISIDEH (752,2) 

COMMON / SOLID / 2 ,R, ARE AE, THSTF 
1 ,THCAP , KPT (752,3) 

REAL*8 ZZ ( 482 , 35 ) , F , EF , B , ELF 1 , ELF2 , ELF3 , FF , HLL 

COMMON / SOLVF.Q / F(482, 1) , EF( 752) 

REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /COST/ KZETA(752) ,KETA(752) ,RHOCP(752) ,BETA(752) 
REAL *8 QRI752) ,HL(752) 

DATA PI/3. 14159265D0/ 

10=6 

TINF = O.DO 


FOR EACH ELEMENT .FORM THERMAL FORCE VECTOR 

WRITE (6, 9999) 

9999 FORMAT (’ FORCE VECTORS ARE ',/) 

DO 110 1=1 ,NE 
EF Cl) = O.DO 

QR(I)=0.0 
110 CONTINUE 

DO 1070 J=1 ,45 
N=2*J 

QR(N)=326 . 0 
1070 CONTINUE 

DO 230 N=1,NE 
ELF(1)=0.D0 
ELF(2)=0 .DO 
ELF(3)=0.D0 
II=NPT(N, 1) 

JJ=NPT(N,2) 

KK=NPT(N , 3) 

Cl = R(KK) -R(JJ) 

CJ = R(II)-R(KK) 

CK = R(JJ)-R(II) 

X(1)=R(II) 

X(2)=R(JJ) 

X(3)=R(KK) 

Y(1)=Z(II) 

Y(2)=Z(JJ) 

Y(3)=Z(KK) 

DO 10 IQ=1 ,2 

IF(ISIDEH(N,IQ).LE.O)GO TO 240 
JQ=ISIDEH (N, IQ) 

KQ=JQ+1 

IF(JQ.EQ. 3) KQ=1 

LG=DSQRT ( (X(KQ) -X(JQ) )**2+ (Y(KQ) -Y( JQ) )**2) 
C WRITE (6. 999 )LG 
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999 FORMAT (2X.F10. 5) 

HL(N)=QR(N)+H(N)*TINF 
IF ( JQ.EQ.1)G0 TO 20 
IF ( JQ. EQ . 2)G0 TO 30 
IF (JQ .EQ.3)G0 TO 40 
GO TO 240 

20 ELF(l) =HL(N)*LG* . 5+ELF ( 1 ) 

ELF(2) =HL (N)*LG* . 5+ELF (2 ) 

GO TO 10 

30 ELF(2) =HL(N)*LG* . 5+ELF (2 ) 

ELF (3 ) =HL (N )*LG* . 5+ELF ( 3 ) 

GO TO 10 

40 ELF(l) =HL(N)*L0*.5+£LF(1) 

ELF ( 3 ) =HL (N ')•• LG* . 5+ELF ( 3 ) 

10 CONTINUE 

C 

240 EF(II)= EF(II)+ELF(1) 

EF(JJ)=EF(JJ) +ELFC2) 

EF(KK)=EF(KK) +ELF(3) 

707 FORMAT (’ ' ,2I8,5F15.5) 

IF(EFCII).EQ.O) GO TO 50 
IF(EF(JJ).EQ.O) GO TO 50 
IFfEF(KK).EQ.O) GO TO 50 
GOTO 230 

C50 WRITE (6, 1999)N,II, JJ,KK,EF(II) ,EF(JJ) ,EF(KK) 

1999 FORMAT ( 1H0 ,4I6,2X,3F10.4) 

230 CONTINUE 

C WRITE (10, 2000) (EF(I), 1=1, NP) 

2000 F0RMAT(2X,F10.4) 

C WRITE (10 ,1000) 

1000 FORMAT (1 HO, 1 LEAVING FORMTF IN TDHEAT ...IB’) 

RETURN 

END 

SUBROUTINE ASSMBL (NP,NBW,ZZ) 

C ASSEMBLE MATRICES FOR RECURRENCE FORMULAS 

REAL *8 Z(482) ,R(4S2) ,AREAE(752) 

REAL*8 DDT , DT,THSTF (482,35), THCAP (482,35) 

COMMON/CONVEN/H (752) , ISIDEH(752, 2) 

C 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 ,THCAP,NPT(752,3) 

REAL*8 ZZ(482, 35) ,F,EF,B ,FF(15 , 1 ) 

COMMON / SOLVEQ / F(482, 1) ,EF(752) 

REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHOCP(752) , BETA(752) 

DO 150 1=1, NP 
C 

C COEFFICIENT MATRIX 

C 

DO 110 J=1,NBW 
ZZ(I.J) = THSTF (I, J) 

110 CONTINUE 

C WRITE(6 , 2000) (ZZ(I,J),J=12,26) 

A '&Q 
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2000 FORMATt’ ' ,6E9.5) 

C 

C THERMAL LOAD VECTOR 

C 

F ( 1 , 1) = EFCI) 

C 

150 CONTINUE 

C WRITE(6, 1313) (F(I , 1) , 1=1 ,NP) 

1313 FORMATC 1 ,G11 .5) 

1001 FORMAT (1H ,13,11011.5) 

C 

C END OF PROGRAM 

C 

RETURN 

END 
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REAL KXX.KXXX 

DIMENSION TITLE ( 10) ,NT1 (II) ,NT2(80) 

DIMENSION NS(3) »T(3) 

DIMENSION I SIDE (2) , IC0N(2) .PHI (3) 

D I MENS ION CFV ( 6 ) , NODNUM ( 6 ) . F I J ( 48 3 , ] 0 0 ) , TI ME ( 1 0 0 ) 

INTEGER NP,IER,N 

REAL*8 Z f 483 ) , R (48 3 ) . AREAE ( 7 1 2 ) , X ( 3 ) , V 1 3 ) , DDT , TIME ,TDT 
2 ,X1 ,X2 ,X3 , Y1 ,Y2 ,Y3 ,COF,THSTF (483,18) .THCAP{483 , 18) .TOLD (483) 
COMMON/CONVEN/H f 752 ) , I SIDEH f 7 5 2 , 2 ) 

COMMON / SOLID / Z ,R , AREAE ,THSTF 
1 , THCAP , NPT (752,3), TOLD 

REAL*8 ZZ(483, 18) ,WK(10000) ,F , EF ,MAX(500) 

COMMON / SOLVEQ / F(483) ,EF(752) ,ZZ 
REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /CONT/ KZETA(752) ,KETA(752) RF0CP(752) ,BETA(752) 

DATA IN/12/,10/6/ 

DATA NT1/343, 347, 351, 355, 359, 363, 367, 371, 375, 379, 383/ 


C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c INPUT OF THE TITLE CARD AND THE CONTROL PARAMETERS 
C<><><><><><><><><><><><><><><><>< ><><><><><><><><><><><><><><><><><> 
C DELTA T: TIME INTERVAL FOR EACH ITERATION 
TDT=500. 0 
DDT=TDT/ 3600 . 

REWIND 12 

C 

RE AD (IN, 3) TITLE 

C WRITE (10,4 ) TITLE 

READ (IN, 2) NP,NE,NBW 


DEFINITION OF THE CONTROL PARAMETERS 

AREAE ( NE ) : ARRAY CONTAINS AREA OF EACH ELEMENTS 
R ( NE*3) : ARRAY CONTAINS THE '.-COORDINATES OF EACH ELEMENT 

Z ( NE*3): ARRAY CONTAINS THE Y-COORDINATES OF EACH ELEMENT 

T ( 3 ) : ARRAY CONTAINS THE TEMP. OF THE ELEMENT’S THREE 

NODES. 

TOLD ( NE*3) : ARRAY COPIES THE T(3) ARRAY FOR EACH ELEMENT 
NE = NUMBER OF ELEMENTS 

NP = NUMBER OF GLOBAL TEMPERATURES 

NBW = BAND WIDTH 
H = CONVECTION COEFFICIENT 
TINF = AMBIENT AIR TEMP. 
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WRITE (10,1102) NP , NE , NBW 

C 

1 FORMAT(4I5 , 6F10 .4) 

2 FORMAT (31 5) 

3 FORMAT (10A4) 

4 FORMAT ( 1H1 , // / IX , 10A4 ) 

1102 FORMAT( IK , 316) 

C 

C 

c 

C OUTPUT OF TITLE AMD DATA HEADINGS 

C 

44 FORMAT (1H1.////1X, 20 A4/ / IX , 5HKXX =,F7 . 1 , 10X, 5HKYY = ,F7.1/1X, 

1 2 1 HC ONVE CTI ON COEFF(l) =,F7. 1,3X,21HC0NVECTI0N C0EFF(2) =,F7.1/ 
2 IX, 15HFLUID TEMP(l) - , F7 . 1 , 3X , 1 5HFLUID TEMP(2) =,F7.1//,1X, 

3 17HNEL NODE NUMBER , 6X , 4HX ( 1 ) , 6X , 4HY ( 1 ) , 6X , 4HX ( 2 ) , 6X , 4HY (2) , 
46X,4HX(3),6X,4HY(3) ) 

C A3SEMBLYING OF THE GLOBAL STIFFNESS MATRIX AND GLOBAL FORCE MATRIX 

C ********* 

C 

C INPUT AND ECHO PRINT OF ELEMENT DATA 
C 

C 

C WRITE (6, 604) 

604 FORMAT (1H1, /////) 

C WRITE (6, 304) 

304 FORMATC T(l) T(2) T(3) 

1 ISIDE(l) ISIDE(2) '/) 

1061 FORMAT (6F10 .3,2110) 

11 FORMATC ' ,3F10.3,2I10) 

C WRITE (6, 1107) 

C WRITE (6, 1108) 

C 

C NE : TOTAL NUMBER OF ELEMENTS 

C 

C . . IN : 12 MAIN LOOP 

T(l)=10 . 

T(2)=10. 

T(3)=10 . 

DO 180 1=1, NE 
H(I)=0 . ODO 
ISIDEH(I , 1)=0 
ISIDEH(I ,2)=0 
180 CONTINUE 

DO 178 JJ=1 ,90 
N=2*JJ 

ISIDEHfN, 1)=2 
ISIDEH(N, 2)=0 
H(N)=5 .DO 
178 CONTINUE 

BIRCH LOOP 






ORIGINAL' PA02 lQ 
do 1001 j=i ,280 POOR QUALITY 

KZETA(J)=.0924 
KETA(J)=. 0836 
RH0CP(J)=10. 94246 

1001 CONTINUE 
C 

C DOUGLAS FIR LOOP 

C 

DO 1002 J=281, 302 
KZETA(J)=.0754 
KETA(J)=.0682 
RH0CP( J)=8 . 5488 

1002 CONTINUE 

DO 1010 J=351 ,638 
KZETA( J)=.0754 
KETA(J)=. 0682 
RH0CP(J)=8 .5488 
1010 CONTINUE 

DO 1015 J=687 , 752 
KZETA(J)=.0754 
KETA(J)=. 0682 
RH0CP(J)=8 .5488 
1015 CONTINUE 
C 

C PAPER LOOP 

C 

DO 1020 J=303 , 350 
KZETA(J)=. 011161 
KETA(J)=. 011161 
RH0CP(J)=8. 1357 
1020 CONTINUE 

DO 1111 J =639,686 
KZETA(J)=.011161 
KETA(J)=. 011161 
RH0CP(J)=8. 1357 
1111 CONTINUE 
C 

C LOOP ON THE REGIONS 

C 

DO 7 KK=1,NE 

IFCKK . GE . 1 . AND . KK . LE . 22 ) BETA (KK)=- . 2618 
IF (KK . GE . 23 . AND . KK . LE . 34) BETA (KK) =~ . 15 707 
IF (KK . GE . 35 - AND . KK . LE . 58) BETA (KK)=0 . 0 
IF(KK.GE . 59 .AND. KK.LE 66) BETA (KK)=. 38394 
IF (KK . GE . 67 . AND . KK . Lb . 74) BETA (KK)= . 17452 
IF (KK . GE . 75 . AND . KK . LE . 82) BETA (KK)=1 . 22164 
IF (KK . GE . 83 . AND . KK . LE . 90 ) BETA (KK)= . 69808 
IF (KK . GE . 9 1 - AND . KK . LE . 98 ) BETA(KK)=2.4433 
IF(KK.GE .99 .AND .KK.LE . 106) BETA (KK)=1 .91972 
IF (KK . GE . 107 . AND . KK . LE . 1 14) BETA (KK)=2 . 9668 
IF(KK.GE. 115. AND. KK.LE. 122) BETA(KK)=2. 7574 
IF (KK . GE . 123 . AND . KK. LE . 146 ) BETA (KK)=3 .14136 
IF (KK . GE . 147 . AND . KK . LE . 158) BETA (KK)=3 . 29843 
IF (KK . GE . 159 . AND . KK . LE . 202) BETA (KK)=3 . 403 14 
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IF(KK.GE . 203 .AND .KK . LE .214) 
IFCKK.GE. 2 15 .AND.KK.LE.230) 
IF(KK . GE . 23 1 . AND . KK . LE . 246) 
IF(KK.GE . 247 . AND. KK . LE . 258 ) 
IF(KK.GE . 259 . AND.KK.LE .326) 
IFCKK.GE . 327 .AND.KK.LE . 350) 
IFCKK.GE. 351. AND.KK.LE. 430) 
IF(KK.GE. 431. AND.KK.LE. 438) 
IF (KK . GE . 43 9 . AND . KK . LE . 446 ) 
IF (KK . GE . 44 7 . AND . KK . LE . 454 ) 
IF(KK.GE. 455. AND.KK.LE. 462) 
IF(KK.GE. 463. AND.KK.LE. 470) 
IF(KK.GE. 471. AND.KK.LE. 478) 
IFfKK.GE. 479. AND.KK.LE. 486) 
IF(KK . GE . 48 7 . AND . KK . LE . 494 ) 
IF(KK .GE . 495 .AND . KK . LE . 502 ) 
IFCKK.GE . 503 . AND.KK.LE .510) 
IF(KK.GE. 5 11. AND.KK.LE. 518) 
IF(KK.GE. 519. AND.KK.LE. 526) 
IF(KK.GE. 527. AND.KK.LE. 534) 
IFCKK.GE.53S. AND.KK.LE. 542) 
IF(KK.GE. 543. AND.KK.LE. 550) 
IFCKK.GE. 551. AND.KK.LE. 558) 
IFCKK.GE. 559. AND.KK.LE. 638) 
IF(KK. GE . 639 . AND . KK. LE . 662 ) 
IFCKK.GE. 663. AND.KK.LE. 708) 
IFCKK.GE. 709. AND.KK.LE. 752) 


BETA (KK)=3. 29843 
BETA(KK)=3. 14136 
BETA(KK)=0.0 
BETA (KK)=-. 15707 
BETA CKK)=-. 2618 
BETA(KK)=-. 15707 
BETA(KK)=0.0 
BETA (KK)=. 38394 
BETA(KK)=. 17452 
BETA (KK)*. 38394 
BETA (KK)*. 17452 
BETA (KK)= 1.22 16 
BETA (KK)=. 69808 
BETA (KK) =1.2216 
BETA (KK)=. 89808 
BETA (KK) =2. 443 3 
BETACKK)=1. 91972 
BETA (KK) =2. 4433 
BETA(KK)=1. 91972 
BETA (KK)=2. 9668 
BETA(KK)=2. 7574 
BETA (KK)=2. 9668 
BETA(KK)=2 .7574 
BETA CKK) =3. 14136 
BETA (KK)=3. 29843 
BETA (KK)=3 . 403 14 
BETA (KK) =0.0 


READ(IN, 1) NEL,NS,XCl),y(l) J XC2) > YC2),X(3),y(3) 

C WRITE(6,121) NEL,NS,X(1) ,Y(1) J X(2),Y(2),X(3),Y(3) ,BETA(KK) 
121 F0RMAT(4I 5 , 7F10 . 4) 

C 

C DIMENSION IN FEET. . . 

X1=X(1) 

X2=X(2) 

X3=X(3) 

Y1=Y(1) 

Y2=Y(2) 

Y3=Y(3) 

C 

X ( 1 )= X(l) /12.0D0 
X ( 2 )= X(2) /12.0D0 
X ( 3 )= X(3) /12.0D0 

Y ( 1 )= Y(l) / 12 . 0D0 

Y C 2 )= Y(2) /12.0D0 

Y ( 3 )= Y(3) /12.0D0 
C 

T1=TC1) 

T2=T(2) 

T3=T(3) 

1107 F0RMAT(1H1 , ///20X 

* T(l) T(2) T(3) 

* / ISIDE(l) ISIDE(2) ') 

1108 FORMAT (' ’ ,20X 


^5 


an oa noon non nano nnooo 
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WRITE (10, 1 1)T, ISIDEH(KK, 1) , ISIDEHCKK, 2) 


TOLD (NS (1))=T(1) 

T0LD(NS(2) )=T(2 J 
TOLDCNSf 3) )=T(3) 

H : CONVECTION COEFFICIENT 

WRITEC6 ,5152)KK,NS (1) ,NS(2) ,NS(3) , TOLD (NS (1 ) ) ,T0LD(NS(2 ) ) , 
1 T0LD(NS(3) ) , H(KK) , (ISIDE(JJ) , JJ=1 ,2) , 

1 XI ,Yl ,X2 , Y2 ,X3 , Y3 ,N0T(KK , 1 ) ,NOT(KK ,6) 

5152 FORMATC > ,415 i 3F10.3 i F5.0,2I2,6F6.4,I5 i A8) 

T(1)=T1 

T(2)=T2 

T(3)=T3 

ECHO THE INPUT 

WRITE (10 ,23) NEL,NS ,X(1) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) 

* , HH,RHO,CP,ISIDE,T 

23 F0RMAT(1H0, 15 , IX, 315 , 2X, 6(1X,F8 . 4) ,2X,2F5 , 0 ,F6 . 3 , 2X,2I3, 

* 2X.3F6.0) 

CALCULATION OF THE CONDUCTION MATRIX 

NPT (NE ,3) : ARRAY STORES THE NODE NUMBERS FOR EACH ELEMENT 
NPT(KK, 1) = NS (I) 

NPT(KK,2) = NS(2) 

NPT(KK,3) = NS (3) 

R(NS(1) ) = X(l) 

R(NS(2)) = X(2) 

R(NS(3) ) = X(3) 

Z(NS(1) ) =Y(1) 

Z(NS(2)) =Y(2) 

Z(NS(3)) =Y(3) 

AREAE (KK)=(X(2)*Y(3) -X(2)*Y ( 1)+X(3)*Y( 1) -X(3)*Y(2) 

1 + X(l)*Y(2)-X(l)*Y(3))/2.0 


WRITE(6,266)KK,X(1) , Y(l) ,X(2) , Y(2) ,X(3) ,Y(3) , AREAE (KK) 
266 F0RMAT( 1H0 , 15 , 6 ( IX , FS . 4 , ) , 2X , D 14 . 7) 

WRITE (6, 3 19 )KK, (T(I) , 1=1,3) 

319 FORMATC' ELEMENT NO. \I8,3F12.4) 

7 CONTINUE 

END OF THE DO LOOP 

WRITE (6 ,264) (I, TOLD (I), 1=1 ,NP) 

264 FORMATC 1H0 , 15 ,E14.5 ,2X, 15 ,E14.5 , 2X, 15 ,E14. 5 ,2X, 15 ,E14.5 

* ,2X, 15 ,E14 . 5 , 2X, 15 ,E14. 5) 

ADD EXTERNALLY APPL. CONC. NODAL FORCES VECTOR TO EF VECTOR 
ID 1=0 
INK=0 

202 READ (IN, 203) NODNUM.CFV 








203 F0RMAT(6I3 , 2X, 6F10 .5) 

ID=0 

DO 204 L=1 , 6 

IF f NODNUM(L).LE.O) GOTO 205 
ID=ID+1 

I NODNL-NODNUM ( L ) 

204 EF {INODNU)=CFVf L)+EF{INODNU) 

GOTO 206 

205 INK-1 
IF(ID.EQ.O) GOTO 216 

206 IFfIDl .EQ. 1)G0T0 222 
WRITE {10, 208) 

208 FORMAT { 1H 1 , ' NODE NUMBER \lOX/ GENERALIZED NODAL FORCES ') 

222 WRITE(I0,207) (NODNUM(L) ,CFV(L) ,L=1,ID) 

207 FORMAT! 1H0 , 6(13 ,E14. 5 , 2X) ) 

IFflNK.EQ. 1)G0T0 216 

I D 1=1 
GOTO 202 
216 CONTINUE 

C 

JJJKNT=0 

C 

DO 3000 ITER=1 , 100 

CALL FORMTK (NE , NP , NBW , TIME , ITER) 

CALL FORMTC (NE.NP.NBW) 

CALL FORMTF (NE,NP,NBW) 

CALL ASSMBL <NP , NBW , DDT, ITER) 

C U r RIT£(6 ,9999) ( (ZZ(I , J) , J=1 ,NBW) ,1=1 ,NP) 

9999 F0RMAT(2X, 7G11 .5) 

C WRITE (6, 9988) (F(I , 1) , 1=1 ,NP) 

9988 FQRMAT(2X,G11 . 5) 

CALL SOLVE (NP, NBW) 

J J JKNT= J J JKNT+ 1 
IF(JJJKNT.GE . 1)G0T0 979 
GOTO 980 

979 CONTINUE 

TIME (ITER)=DDT*ITER 
DO 2001 I = 1,483 
FIJ(I , ITER) = F(I) 

2001 CONTINUE 

C WRITE ( 10 , 666 )TIME , ITER , TDT 

C66 FORMATC// , 1H1 , ' TEMPERATURE DISTRIBUTION AT THE TIME =' , 

C * E15.6, ’ HOUR ’ ,/' ITERATION NO. =’,110, *#*, 

C * * *&* TDHEAT.2 DDT= ' , F 10 . 6 , / ) 

C IF(ITER.EQ.22) 

C 1 WRITE (7 ,777)(I,F(I),I=1,NP) 

C77 FORMATC 1H0 , 15 ,F14. 5 , 2X, 15 ,F14.5 ,2X, 15 ,F14 . 5 ,2X, 15 ,F14. 5 , 

C *2X, 15 ,F14 . 5 , 2X, 15 ,F14.5) 

C77 FORMAT { 14, F 10. 4) 

980 CONTINUE 
MAX(1)=0.0 
STT=0 . 0 

DO 1515 1=1, NP 

A 'll 




ORKSM. PAvS ‘3 
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MAX C 1+1 )=DABS (F ( I ) -TOLD (I ) ) 

CF WRITE (6,121 1 )MAX ( 1+1 ) 

1211 FORMAT ( 2X , F IS . 6 ) 

IF (MAX (1+1 ) . LT . MAX ( I ) ) GOTO 1515 
IF (MAX ( I + 1 ) . GT . MAX ( I ) ) FAXM=MAX ( I +1 ) 

IF(FAXM.GT. STT)GTST=FAXM 
STT=GTST 
1515 CONTINUE 

C WRITE(6, 1222)GTST 

CF GTST IS THE STOPPING CRITERION. ITERATION STOPS WHEN DIFF. BETW 
CF T(NEW) AND T(OLD) IS LESS THAN 0.3 DEGREE 
1222 FORMAT (/ ,2X,F10 .5) 

IF (GTST. LT. 0.3) GOTO 1200 
DO 2000 1=1, NP 
TOLD(I)=F(I) 

2000 CONTINUE 

3000 CONTINUE 

1200 DO 3011 IJP = 1,483 

WRITE (7, 666) IJP 

666 F0RMAT(//, ' TEMPRATUE DISTRIBUTION AT POINT ',15,/) 

C WRITE (7 , 777) (I , TIME (I) ,FI J(I JP , I) , I =1,65) 

777 FORMAT (I5,2F10.4) 

3011 CONTINUE 

C 

C 

C200 STOP 


END 

SUBROUTINE FORMTK (NE,NP,NBW, TIME, ITER) 

Q****f-*V:V:***************V:*V;**VfVf***V-*VfVt**-,V****Vt***V.-***‘.’;‘.V****‘,T-,'rVc**-.'j**'.V-'.--.V 


c. 

c 


FORM BANDED THERMAL CONDUCTANCE MATRIX 


IMPLICIT REAL*4(A-H,0-Z) 

REAL*8 KZETA , KETA , RHOCP , BETA 

REAL*8 B I , B J , BK , Cl , C J , CK , COF , CONI , C0N2 , CON3 

REAL*8 RBAR , CON , LG , RCTL , ELKT (3 , 3) , THSTF (483 , 18 ) 

REAL*8 A00 ,A11 ,A22 ,B00 , Bll ,B22 ,C00 ,THCAP(483, 18) , TOLD (483) 
REAL*8 Z(483) ,R(483) ,AREAE(752) ,X(3) ,Y(3) , DDT, TIME ,ATIM,P 
COMMON/CONVEN/H ( 75 2 ) , I S IDEH (752,2) 


C 

C 

C 

C 


COMMON / SOLID / Z,R, ARE AE, THSTF 
1 ,THCAP,NPT(752,3),TOLD 
REAL*8 ZZ(483, 18) ,F, B,FF,EF 
COMMON / SOLVEQ / F(483) ,EF(752) ,ZZ 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RH0CP(752) , BETA(752) 


INITIALIZE GLOBAL THERMAL STIFFNESS MATRIX 

IN=5 


10=6 


DO 120 J=1 ,NBW 
DO 110 1=1, NP 
THSTF ( I, J) = 0.D0 
110 CONTINUE 


120 CONTINUE 


OF POO!? QukU» 


c 

c 

130 NBAND = 0 

DO 899 J=1 ,NE 

C WRITE(6 , 788)KZETA(J) ,KETA( J) ,BETA(J) 

788 FORMAT ( 3X , 3F20 . 9 ) 

899 CONTINUE 
C 

C WRITE (6 ,7776) 

7776 FORMAT (' CONDUCTION MATRIX ELEMENTS ARE ' ,/) 

DO 230 N=1 ,NE 

C<><> <><><><> <><>«:><><><><><><><><> 

C 

II = NPT(N.l) 

JJ = NPT(N, 2) 

KK = KPTCN.3) 

BI = Z( JJ)-Z(KK) 

BJ = Z (KK) -Z (II ) 

BK = Z(II)-Z(JJ) 

Cl = R(KK)-R(JJ) 

CJ = R(II)-R(KK) 

CK = R(JJ) -R(II) 

X(l) =RCII) 

X(2) =R( JJ) 

X(3) =R(KK) 

YC1) =Z(II) 

Y(2) =Z(JJ) 

Y(3) =Z(KK) 

C WRITE(6,555)X(1) ,X(2) ,X(3) ,Y(1) ,Y(2) ,Y(3) 

555 FORMAT ( 2X , 6F12 . 6 ) 

ALPHA=BETA(N) 

CON 1= ( KZETA (N ) * CCOS (ALPHA) )**2+KETA (N)* (SIN (ALPHA) )**2) 
1/(4.*AREAE(N)) 

CON2=(KZETA(N) -KETA(N) )*SIN ( 2 . --ALPHA)/ (4 . * AREAS (N) ) 

C0N3= ( KZETA ( N ) * ( S IN (ALPHA ) ) **2+KETA ( N ) * ( COS (ALPHA ) ) **2 ) 

1/(4 .*AREAE (N) ) 

C 

C WRITE (6 , 799 9 ) CONI , CON2 , COM3 , AREAE (N ) 

7999 F0RMAT(2X,3F10.3,2X,F10.7) 

ELKT (1,1) = (B 1**2 )*C0N1+ (CI**2 ) *C0N3 
ELKT(1 ,2) = ( B I*B J ) *C0N 1+ ( C I*C J ) *C0N3 
ELKT(1,3) = ( B I *BK ) *CON 1+ ( C I *CK ) *C ON3 
ELKT (2,2) = ( B J-'"- 2 ) *C0N1+ ( C J**2 ) *C 0N3 
ELKT(2,3) = (BJ*BK)*C0N1+(CJ*CK)*C0N3 
ELKT(3 ,3) = ( BK**2 ) *CON 1+ (CK**2 ) *CON3 
DO 10 IQ=1,2 

C WRITE (6, 7805 )N,( ISIDEH(N,J1) , Jl=l,2) 

7805 FORMAT (’ SIDE EXPOSED ’,3110) 

IF ( ISIDEH(N, IQ) . LE . 0 ) GO TO 240 
JQ = ISIDEH (N, IQ) 

C WRITE (10 ,12) JQ,N 

12 FORM AT (1 HO , ' INSIDE' ,13 ,3X, ' OF ELEMENT' ,15) 

KQ = JQ +1 


no no n n oononoono 
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IF ( JQ.EQ . 3) KQ=1 

LG = DSQRT( (X(KQ) -X(JQ) )**2+(Y (KQ) -Y ( JQ) )**2) 

HL =H(N)*LG/6. 

IFfJQ.EQ.l) GO TO 20 
IFfJQ.EQ.2) GO TO 30 
IFCJQ.EQ.3) GO TO 40 
GO TO 10 

20 ELKT(1 , 1 )= ELKT( 1 , 1 )+HL * 2. 

ELKT(2,2)= ELKT(2,2)+HL * 2 . 

ELKTfl ,2)= ELKTf 1,2)+HL 
GO TO 10 

V 

30 ELKT(2,2)= ELKT(2,2) J -HL*2. 

ELKT(3 , 3)= ELKT(3.3 )+HL*2. 

ELKT(2,3)= ELKT(2 ,3)+HL 
GO TO 10 

40 ELKT(1 , 1)= ELKT f 1 , 1 j +HL * 2. 

ELKT(3,3)= ELKTf3 ,3)+HL * 2. 

ELKTfl ,3)= ELKTfl ,3)+HL 
10 CONTINUE 

WRITE fIO,2000)N, II jJJjKK,X( 1) , Yfl) ,Xf2) ,Y(2) ,Xf3) ,Y(3) ,CONDTY (N) , 
lAREAEfN) , RBAR , CON , ELKTf 1,1) , ELKT f 1 , 2 ) , ELKT (1,3) ,ELKT(2,2) , 
2ELKTf2,3) ,ELKTf3 ,3) 

2000 FORMATf 1HO , 415 , 8G1 1 . 4 , /2 IX , 8G1 1 . 4) 

Oj 

STORE IN GLOBAL MATRIX c{ 

01 

CONTINUE 

WRITEf I0,266)N,Xf 1) ,Yf 1) ,Xf2) ,Yf2) ,Xf3) ,Yf3) ,AREAE(N) , NB AND 
266 FORMAT f 1H0 , 1 5 , 6 f IX , F8 . 4) ,2X,D14. 7 ,2X, 15) 

WRITEf 6, 277)N,(f ELKTf I , J) , J=1 , 3) , 1=1 , 3) 

277 FORMATf’ \I5/ 3f3F25.6/)) 

DO 220 LL=1,3 
DO 210 MM=1 , 3 
IFfMM.LT. LL) GOTO 210 
I=NPTfN,LL) 

J=IABS fNPTfN,MM) -I)+l 

IFfNPTfN,M>l) .LT.NPTfN,LL)) I=NPTfN,MM) 

THSTF f I , J) =THSTF (I, J)+ELKT(LL, MM ) 

WRITE f6 , 7543 ) I , J , 1 1 , JH , NPT f N , J) , NPT f N , I ) 

7543 FORMATf ',618, ’******') 

IF (J.GT.NBAN'D) NBAND=J 
210 CONTINUE 
220 CONTINUE 
230 CONTINUE 

WRITEf6 ,5050) f fTHSTFf I , J) , J=1 ,NBW) , 1=15,20) 

5050 FORMATf 7G11. 5) 

WRITE (6,200) NBAND 
200 FORMATf lHOf NBAND= ’,15) 

END OF PROGRAM 


A 
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RETURN 

END 


FORM BANDED CAPACITANCE MATRIX FOR A DEFORMING F.E. MESH 

SUBROUTINE FORMTC (NE.NP.NBW) 

REAL*8 Z(4S3) ,R(483) ,AREAE(752) ,X(3) , Y ( 3 ) , DDT 
REAL * 8 CON , ELCAP f 3,3) ,THCAP(4S3 ,18) , TOLD (48 3) ,THSTF(483 , 18) 
1 ,THSTF1(747, 17) 

COMMON/ CONVEN/H ( 75 2 ) , S IDEH ( 7 5 2 , 2 ) 

COMMON / SOLID / Z , R , ARE AE , TH STF 
1 ,THCAP,NPT(752 ,3) ,TOLD 
REAL--8 ZZ (483 , 18 ) , F , EF , B , COF , FF , KZETA , KETA 

COMMON. / SOLVEQ / F(483) ,EF(752) , ZZ 
DIMENSION FF ( IS , 1 ) 

COMMON /CONT/ KZETA(752) ,KETA(7S2) ,RHOCP(752) ,BETA(752) 


REAL * 8 R1 , R2 , R3 , R4 , R5 , R6 

INITIALIZE GLOBAL THERMAL CAPACITANCE MATRIX 


10=6 

J9=0 

DO 120 J=1,NBW 
DO 110 1=1, NP 
THCAP(I,J) = 0 .DO 
110 CONTINUE 
120 CONTINUE 


WRITE (6, 8880) 

80 FORMAT ( ' CAPACITANCE MATRIX ELEMENTS ARE ',/) 

200 DO 230 N=1 ,NE 

FOR EACH ELEMENT, FOR THERMAL CAPACITCANCE MATRIX 

II = NPT(N, 1) 

JJ = NPT(N,2) 

KK = NPT(N,3) 

CON=RHOCP (N ) *AREAE (N ) / 1 2 . 

ELCAP ( 1 , 1 ) =2 . *CON 
ELCAP ( 1 , 2 ) =1 . *CON 
ELCAP (1,3)=1.*C0N 
ELCAP(2, 2)=2.*C0N 
ELCAP ( 2 , 3 ) =1 .*C0N 
ELCAP ( 3 , 3 )=2 . *CON 


W t RITE(IO,2000)N,II,JJ ) KK,R(II;,Z(II) ,R(JJ) ,Z(JJ) ,R(KK) ,Z(KK) , 


non o nonon non non oo no no 
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1 RHOCP(N) ,AREAE(N) ,C0N,ELCAP(1 , 1) ,ELCAP(1 ,2) ,ELCAP( 1 ,3 ) , 

2 ELCAP(2,2) ,ELCAP(2,3) ,ELCAP(3,3) 

2000 F0RMAT(1H0,4I5 ,9G11 .4/21X.6G11 .4) 

STORE IN GLOBAL MATRIX 

DO 220 1=1,3 
DO 210 J=1 , 3 
IF(J.LT.I) GOTO 210 
II=NPT(N,I) 

JJ=IABS(NPT(N,J)-11)+1 
IF (NPT(N,J).LT.NPT(K,I)) II=NPT(N,J) 

THCAP(II,JJ) = THCAP(II,JJ)+ELCAP(I,J) 

210 CONTINUE 

220 CONTINUE 

WRITE (6, 277 )N, (( ELCAP(I , J) ,J=1 ,3) , 1=1 ,3 ) 

277 F0RMATC ' ,15/ 3(3F25.6/)) 

230 CONTINUE 

WRITE ( 10 , 4000 ) (THCAP( I , J) , J=1 , NBW) 

4000 FORMAT ( ] HO , ' THC AP ' , 10G 1 1 . 4/ 6X , 1 0 G 1 1 . 4 / 6X , 1 OG U . 4 / 6X , 1 0G 1 1 . 4 ) 
WRITE(I0, 1000) 

1000 F0RMATC1H0, 'LEAVING F0RMTC IN TDHEAT ..IB ') 

264 FORMAT ( 1H0 ,15 } E 14.5 ,2X, 15 ,E14.5 ,2X,I5 ,E14.5 ,2X, 15 ,E14.5 

* ,2X,I5 ,E14 .5 ,2X,I5 ,E14.5) 

END OF PROGRAM 

RETURN 

END 

SUBROUTINE F0RMTF ( NE , NP , NBW) 

FORM THERMAL FORCE VECTOR 


REAL* 8 ELF (3) ,LG,RCTL,C0F,THCAPC4S3 , 18) ,THSTF(483 , 18) 
REAL*S Z (483) ,R(483) ,AREAE(752) ,X(3) ,Y(3) , DDT, TOLD (483) 
COMMON/ CON VEN/H ( 7 5 2 ) , I S IDEH ( 7 5 2 , 2 ) 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 ,THCAP,NPT(752 ,3) ,T0LD 

REAL*8 ZZ (483,18) ,F,EF,B,ELF1 ,ELF2,ELF3 ,FF,HLL 

COMMON / SOLVEQ / F(483) ,EF(752) , ZZ 
REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHOCP(752) . BETA(752) 
REAL *8 QR(752) ,HL(752) 

DATA W / .666666666600/ 

10=6 

TINF =10. DO 


FOR EACH ELEMENT .FORM THERMAL FORCE VECTOR 


ORIGINAL £>&&£ q 
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C WRITE(6 , 9999 ) 

9999 FORMAT (' FORCE VECTORS ARE ’ ,/) 

DO 110 1=1, NE 
EF(I) = O.DO 
QR(I)=0 . 0 
110 CONTINUE 

DO 1070 J=1 ,45 
N=2*J 

QR(N)=326 . 0 
1070 CONTINUE 

DO 230 N=1 ,NE 
ELF (1 )=0 .DO 
ELF(2)=0 .DO 
ELF(3)=0 .DO 
II=NPT(N,1) 

JJ=NPT(N,2) 

KK=NPT(N ,3 ) 

Cl = R(KK) -R(JJ) 

CJ = R(Ilj-RCKK) 

CK = R(JJ)-R(II) 

BI=Z(JJ) -Z(KK) 

BJ=Z(KK) -Z(II) 

BK=Z(II) -Z(JJ) 

X(1)=R(II) 

X(2)=R(JJ) 

X(3)=R(KK) 

Y(i)=zcn) 

Y(2)=ZCJJ) 

Y(3)=Z(KK) 

ALPHA=BETA(N) 

C0N2= (KZETA(N) -KETA (N) )*SIN ( 2. -'ALPHA) / (4 . *AREAE (N) ) 

ELF 1=C0N2* (B I*C I -TOLD (II) +B I*C J*TOLD ( J J ) +B I*CK*TOLD ( KK ) ) 

ELF 2=C0N 2* ( B J*C I*TOLD ( 1 1 ) +B J*C J-TOLD (JJ)+BJ*CK*TOLD (KK) ) 

ELF3=CON2* (BK*CI*TOLD ( 1 1 )+BK*C J*TOLD (JJ)+BK*CK*TOLD (KK) ) 

ELF ( 1 )=-W*ELFl 

ELF(2)=-W*ELF2 

ELF ( 3 )=-W*ELF3 

DO 10 IQ=1 } 2 

IF(ISIDEH(N, IQ) .LE.O)GO TO 240 
JQ=ISIDEH (N, IQ) 

KQ=JQ+1 

IF(JQ.EQ.3) KQ=1 

LG=DSQRT( (X(KQ) -X( JQ) )**2+ (Y(KQ)-Y(JQ) )**2) 

C WRITE(6,999)LG 

999 FORMAT (2X,F10. 5) 

HL (N ) =QR (N ) +H (N)*TINF 
IF ( JQ.EQ. 1)G0 TO 20 
IF ( JQ.EQ. 2)G0 TO 30 
IF (JQ .EQ.3)G0 TO 40 
GO TO 240 

20 ELF ( 1 ) =HL(N)*LG* . 5+ELF ( 1 ) 

ELF(2) =HL(N)*LG* . 5+ELF (2 ) 

GO TO 10 

30 ELF(2) =HL ( N ) *LG* . 5+E LF ( 2 ) 


noon 
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40 

10 

C 

240 


707 


C50 

1999 

230 

C 

2000 

C 

1000 


C 


C 


C 

c 


c 


c 

1000 

c 

6567 

3000 

4000 

C 


ELF(3) =HL (N ) *LG* . 5+ELF ( 3 ) 

GO TO 10 

ELF(l) =HL(N )*LG* . 5+ELF ( 1 ) 

ELFI3) =HL (N) *LG* . 5+ELF f 3 ) 

CONTINUE 

EF ( 1 1 ) = EF(II )+ELF(l) 

EF(JJ)=EF(JJ) +ELF (2) 

EF(KK)=EF(KK) +£LF(3) 

FORMAT ( ' ' .2I8.5F15.5) 

IF(EFCII).EQ.O) GO TO 50 
IF(EF(JJ) .EQ.O) GO TO 50 
IFCEF(KK).EQ.O) GO TO 50 
GOTO 230 

WRITE (6 , 1999 )N, 1 1 . JJ ,KK ,EF( II) ,EF( JJj ,EF(KK) 

FORMAT( IHO ,416 , 2X, 3F10 . 4) 

CONTINUE 

WRITE (10, 2000) (EF(I), 1=1. NP) 

FORMAT ( 2X , F 1 0 . 4 ) 

WRITE (10 ,1000) 

FORMAT (IHO,' LEAVING FORMTF IN TDHEAT ...IB 1 ) 

RETURN 

END 

SUBROUTINE ASSMBL (NP ,NBW ,DDT, ITER) 

ASSEMBLE MATRICES FOR RECURRENCE FORMULAS 

REAL*8 Z (483) ,R(483) ,AREAE(752) ,THSTF(483 , 18) ,THCAP(483 , 18) 
. ,THSTF1(747 , 17) 

REAL*8 ZZ(483, 18) ,F,EF,B,ZO(747, 17) ,FF(15 , 1) ,T0LD(483) 
COMMON/CONVEN/H(752) , ISIDEH(752 ,2) 

COMMON / SOLID / Z , R , AREAE , THSTF 
,THCAP,NPT(752,3) ,TOLD 
COMMON / SOLVE Q / F(483) ,EF(752) ,ZZ 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHOCP(752) ,BETA(752) 

WEIGHTING FACTORS FOR GALERKIN ANALYSIS 

DATA W1 ,W2 ,W3 , W4/ 0.333333333333333D0 , 

0 . 50000000000000000 , 

! -0.1666666666666667D0 , 

i 0 . 500000000000000DQ / 

10=6 

DT = DDT 

WRITE (10,1000) W1 I V2,W3,V4,DT 

FORMAT (IHO, 'ASSEMBLE IN TDHEAT 1' // ,5E15 .8/' ZZ 

. COL 1 TO 6 '/) 

FORMATC 1 , 15 , 6F12. 6) 

FORMAT (IHO, 'THSTF* , 10G11 ,4/ 6X, 10G11 -4/6X, 10G11 . 4/6X, 10G11 .4/) 
FORMAT ( IHO , ’ THCAP ' , 1 OG 1 1 . 4/ 6X , 1 0G1 1 . 4/6N , 1 OG 11 . 4/ 6X , 10G1 1 . 4/ ) 
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DO 150 1=1, NP 
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C COEFFICIENT MATRIX 

C 

DO 110 J=1,NBW 

ZZ(I,J) = Wl* (THSTF (I , J) )+W2*(THCAP(I , J) )/DT 
110 CONTINUE 

C WRITE (10, 2000) (ZZ(I , J) , J=1 ,NBW) 

2000 FORMATC ' , 8E15 .6) 

THERMAL LOAD VECTOR 

F( I ) = EF(I)/2.D0 

IF (I.EQ.l) GO TO 130 
JST = 1+1 -NBW 
IF (JST.LT.l) JST=1 
JEND = 1-1 

WRITE (6,1000) I. JST, JEND 

DO 120 J= JST, JEND 
II = I+l-J 

C! 

F ( I ) =F(I)+ (W3* (THSTF ( J , 1 1 ) )+W4* (THCAP ( J , 1 1 ) ) 

1 /DT)* (TOLD(J) ) 

WRITE (6,1000) I,J,II,THSTF(J,II),THCAP(J,II),TOLD(J),F(I) 
120 CONTINUE 

130 JST = I 

JEND = I-l+NBW 
IF (JEND.GT.NP) JEND=NP 
WRITE (6,1000) JST, JEND 
DO 140 J=JST, JEND 
JJ = J-I+l 

F(I) =F ( I )+(W3* (THSTF ( I , JJ) )+W4*(THCAP(I , JJ) ) 

1 /DT)* (TOLD(J)) 

WRITE (6,1000) I , J, JJ,THSTF(I , JJ) , THCAP (I , JJ) ,T0LD( J) ,F(I ) 
140 CONTINUE 

WRITE (6,1000) I, JST,JEND,F(I) 

WRITE (6,1001) I,(ZZ(I,J),J=1,NBW),F(I),T0LD(I) 

150 CONTINUE 

1001 FORMAT (1H ,I3,6G11.5) 

END OF PROGRAM 

RETURN 

END 

SUBROUTINE SOLVE (NP ,NBW) 

C 

C SOLUTION OF EQUATIONS BY GAUSS ELIMINATION 

C NP=NO. OF NODAL POINTS N'B W=B ANDW I DTH 

C =NO. OF UNKNOWNS 

C =N0. OF EQUATIONS 

C 

REAL*8 ZZ ( 48 3 , 1 8 ) , F , EF , B , P , EXTRA 
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COMMON/SOLVEQ/F (483) ,EF(752) ,ZZ 
10=6 

C WRITE (6, 1002) 

1002 FORMAT ( ' IN SOLVE ZZfl.J) COL 1 TO 4 '/) 

100 DO 130 1=1, NP 

C WRITE! 10 , 1001 ) I , (ZZ(1 , J) , J=1 ,4) ,F(I ) 

1001 FORMAT! ' ' ,I5,3E15.5) 

REDUCE THE STIFFNESS MATRIX 

DO 120 J=2,\Blv 
II=I+J-1 

IF(ZZ(I , J) .EQ.0.0) GOTO 120 

p= zza,j)/;:zn,n 

JJ=0 

DO 110 K=J,\BW 
JJ=JJ+1 

IFCZZf I ,K) AE.0.0D0) ZZ(II ,JJ)=ZZ(II ,JJ)-P*ZZ(I ,K) 

IF(DABS(ZZ( 1 1 , JJ)) .LT. 1 .E-20)ZZ(II ,JJ)=0.0D0 
110 CONTINUE 

ZZ(I,J) =P 

REDUCE LOAD VECTOR 

F(II) = F(II)-P*F(I) 

120 CONTINUE 

F(I) = F(I)/ZZ(I,1) 

WRITE CIO, 1000)1 , ZZ(I , J) , J=1 ,NBW) ,F(I) 

1000 FORMAT (1 HO, ' SOLVE '.I 3/4 C5E 15. 8/)) 

130 CONTINUE 

BACK SUBSTITUTION 

200 N=NP 

210 N=N-1 

IF(N.EQ.O) GO TO 300 
II=N 

DO 220 J=2.NBW 
11=11+1 

IF ( ZZ(N, J) .NE. 0 . 0) F(N)=F(N) -ZZ(N, J)*F(II) 

220 CONTINUE 

GO TO 210 

SOLUTION VECTOR STORED IN F ARRAY 

300 RETURN 

END 
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C 

C THE GRID PROGRAM IS MODIFIED TO GENERATE RECTANGULAR ELEMENTS. 

C THE CHANGES ARE AS FOLLOW: 

C 

CH THE DIMENSION OF N00(2500,6) IS CHANGED TO N00(2500,7) 

CH THE DIMENSION OF NR (3) IS CHANGED TO NR (4) 

CH THE DIMENSION OF IS CHANGED TO LB(3) TO LB(6) 

C 

C DIVISION INTO RECTANGULAR ELEMENTS 
C 

K=1 

DO 17 1=1 jNROWS 
DO 17 J=1 , NCQL 
XE(K)=XC(I , J) 

YE(K)=YC(I,J) 

NE(K)=NN(I , J) 

17 K=K+1 
L=NROWS-l 
DO 21 1=1, L 
DO 21 J=2,NCOL 

CH DIAG1=SQRT( (XC (I , J) -XC ( 1+1 , J- 1 ) )**2+(YC (I , J) -YC (1+1 , J- 1 ) )**2) 

CH DIAG2=SQRT( (XC (1+1 , J) -XC ( I , J- 1 ) )**2+(YC ( 1+1 , J) -YC ( I , J- 1) )**2) 

NR( 1 )=NCOL*I+J- 1 
NR(2)=NCOL*I+J 
NR(3)=NC0L*(I-1)+J 
NR (4 )=NCOL* ( I - 1 )+J- 1 
CH DO 21 I J=1 , 2 
NEL=NEL+1 

CH IF((DIAG1/DIAG2) .GT. 1.02) GO TO 18 
J1=NR(1) 

J2=NR(2) 

J3=NR(3) 

J4 = NR (4) 

CH GO TO 19 

CH 18 J1=NR(IJ) 

CH J2=NR(IJ+1) 

CH J3=NR(4) 

CH LB (3), LB (S) AND LB (6) ARE ADDED 

LB ( 1 )=IABS (NE ( J1 ) -NE ( J2) )+l 
LB (2)=IABS (NE ( J1 ) -NE ( J3 ) )+l 
LB (3)=IABS (NE ( J1 ) -NE ( J4 ) )+l 
LB (4)=IABS (NE ( J2) -NE ( J3 ) )+l 
LB(S)=IABS (NE (J2) -NE(J4) )+l 
LB (6)=IABS (NE ( J3 ) -NE ( J4) )+l 
DO 20 IK=1 , 6 

IF (LB (IK) . LE . NBW) GO TO 20 
NBW=LB(IK) 

NELBW=NEL 
20 CONTINUE 

IF( KK . NE . 46 ) GO TO 888 
C WRITE(6,666) J3 , NE(J3) 

666 FORMAT (///// ' J3 = ',19,' NE(J3) = ' , 19 ) 

888 CONTINUE 

WRITE(IO, 113) KK,NEL,NE(J1) ,NE(J2) ,NE(J3) ,NE(J4) ,XE(J1) ,YE(J1), 
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1 XE ( J2 ) , YE ( J2) , XB ( J3 ) , YE ( J3 ) , XE ( J4) , YE ( J4) 

IF(IPCH.EQ.O) GO TO 21 

WRITEUP, 114) NEL,NE(J1) ,NE(J2) ,NE(J3) ,NE(J4) ,XE(J1) ,YE(J1) ,XE(J2) 
1 , YE(J2) ,XE(J3) ,YE(J3) ,XE(J4) , YE(J4) 

CH NOO(NEL, 6) IS ADDED 
NOO(NEL, 1)=KK 
N00{NEL,2)=NEL 
N00(NEL,3)=NE(J1) 

N00(NEL,4)=NE(J2) 

NOO (NEL , 5 ) =NE ( J3 ) 

NOO ( NE L , 6 ) =NE ( J4 ) 

NOO(NEL,7)=IRR(KK) 

21 CONTINUE 

22 CONTINUE 
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CH THE DIMENSION OF NOO, X, I, OR Y IS INCREASED BY 1 
C 

COMMON/ AAA/ NOO ( 25 00 , 7 ) , NTT ( 8 0 ) 

COMMON/BBB/ X(4,2000) ,Y(4,2000) ,NBO(100G) ,TITLE(10) 
COMMON/CCC/ NODES ,LMENTS, JT( 12000) ,MEMJT(24000) , JMEM(3000) , 
$JNT(3000) , ID IFF, M INMAX 
DATA IN/ II/, 10/ 6/, IP/ 12/ 

C 

READ (IN, 100)TITLE 
WRITE(IP, 100)TITLE 
100 FORMAT(10A4) 

122 F0RMAT(I3 .5F10.5) 

J =1 
JJ=0 

CH I IS CHANGED FROM 3 TO 4 

150 READ (IN, 17) (JT(3000-' f (I- 1)+J) , 1=1,4) , 1(1, J) ,Y(I,J) ,1=1,4) 

CH JT(9000+J)=0 

C WRITE (10, 17) (JT(3000*(I -1 )+J) , 1=1,4) , ("(I , J) ,Y(I , J) , 1=1,4) 
NODES=MAXO ( JT ( J ) , JT(3000+J) , JT(6000+J) , JT(9000+J) , JJ) 
JJ=NODES 

IF(JT(3000+J).EQ.O)GO TO 152 
J=J+1 
GO TO 150 
152 LMENTS=J-1 
NODES=JJ 

12 FORMAT (/ / , 5X, 15HNCJMBER OF NODES, 14, 15X, 

1' NUMBER OF ELEMENTS ’,14,//) 

WRITE (10, 105) TITLE 
105 FORMAT(//,1X,20A4,//) 

C WRITE (IQ , 12)N0DES , LMENTS 

C WRITE (10, 13) 

C3 FORMAT(// , 3X, 18HNEL JT1 JT2 JT3 , 1 0X, -HX(1) , 8X, 

C $4HY(1) ,8X,4HX(2) ,8X,4HY(2) ,8X,4HX(3j ,8X,4HY(3)) 

17 F0RMAT(4X,4I4,8F10 .5) 

DO 10 J=l, LMENTS 

WRITE(IO, 11)J, (JT(3000*(I-1)+J) ,1=1,4) , (X(I ,J) ,Y(I , J) ,1=1,4) 
11 F0RMAT(1X,5I5 ,3X,8F12 .4) 

10 CONTINUE 

READ (IN, 300) IBO 
300 FORMAT (5 OX, 14, 46X) 

READ(IN,301) (NBO(I) , 1=1, IBO) 

C WRITE (10, 302) 

DO 333 1=1, LMENTS 

READ (IN, 444) (NOO(I , J) , J=1 ,7) 

333 CONTINUE 

C WRITE(IO,303)(NBO(I),I=1,IBO) 

302 FORMAT(//,’ *** BOUNDARY NUMBERS ***’ ) 

CALL SETUP 
NTBAN=ID IFF+ 1 
WRITE (10, 202)KTBAN 

202 FORMAT (/ / , 2X , 25HTHE ORIGINAL BANDWIDTH IS, 14,//) 

CALL OPTNUM 

IF ( IDIFF . LE . MINMAX) GO TO 115 
MIBAN=MINMAX+1 
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WRITE(IO, 198) 

198 FORMATC 1H1/// , IX , 5 (3H0LD , 3X , 3HNEW, 7X) , / , IX , 

$5 (4HN0DE , 2X , 4HN0DE , 6X) ) 

WRITE ( 7 , 200 ) ( J , JNT ( J) , J=1 , NODES ) 

200 FORMAT(5 (IS , IX, 15 ,5X) ) 

C WRITE (10, 13) 

DO 180 J=1 , LMENTS 
DO 555 11=1,4 
I I 1=11+2 

555 N00(J,III)=JNT(JT(3000*(1I-1)+J)) 

C 

180 CONTINUE 

C WRITE (10, 11) J, (JNT(JT(3000*(I-1)+J) ) , 1=1,4) , 

C $(X(I , J) ,Y(I,J) ,1=1,4) 

WRITE (10, 201 )MIBAN 

201 FORMATC // , 2X , ' THE NEW BANDWIDTH IS ' , I 3 , / ) 

C WRITE ( I P , 3 10 ) NODES , LMENTS , MI BAN 

C DO 181 J=l, LMENTS 

C81 WRITEUP, 18) J, (JNT(JT(3000*(I-1)+J)) ,1=1,4) , (X(I,J) , 

C $Y(I,J) ,1=1,4) 

08 F0RMAT(5I5 ,8F10 .4) 

C WRITE (10, 302) 

C WRITE (IP , 300) IBO 

C WRITE ( IP , 30 1 ) (JNT (NBO (I)), 1=1, IBO) 

C WRITE (10,303) (JNT (NBO ( I ) ) , 1=1 , IBO) 

303 FORMAT (15 14) 

GO TO 117 
115 CONTINUE 

WRITE(IO, 101) 

117 CONTINUE 

101 FORMATC///, 2X, ’THE ORIGINAL BANDWIDTH IS MININMUM ’ ) 

WRITE (IP , 3 10 )NODES , LMENTS , NTBAN 
310 FORMAT (315) 

C DO 190 J=l, LMENTS 

C190 WRITE(IP, 18) J, (JT(3000*(I-1)+J) ,1=1,4) , 

C 1(X(I,J),Y(I,J),I=1,4) 

C WRITE (IP, 300) IBO 

C WRITE (IP, 301) (NBO (I) ,1=1, IBO) 

301 FORMAT (2015) 

C Cl 

444 FORMAT (615, A4) 

LLL = NOO (1,7) 

C 

LOL=0 

MMM=0 

NNN=1 

DO 666 11=1, LMENTS 

IF ( LLL. EQ. NOO (II, 7) )G0 TO 1888 

NNN=NNN-1 

LOL=LOL+NNN 

MMM=MMM+1 

NTT(MMM)=LOL 

C WRITE (10, 888) NNN,NTT(MMM) 

LLL=N00(II ,7) 


NNN=1 

1888 CONTINUE 

WRITE (10 , 777)NNN, (NOD (II , J) , J=l,7) , (X(I,II),Y(I,II), 1=1,4) 
C WRITE (7 , 1012)N00(II ,2) ,NOO(II ,6} 

WRITE C 7 ,199) (NOO(IIjJ) ,J=2,6), (X(I,II),Y(I,II) ,1=1,4) 

IF (IT. NE . LMENTS ) GO TO 1889 

LOL=LOL+NNN 

MMM=MMM+1 

NTTCMMM)=LOL 

C WR ITE ( 1 0 , 8 8 8 ) NNN , NTT ( MMM) 

1889 NNN=NNN+1 
666 CONTINUE 

C 

C1012 F0RMAT(I5,A5 ,IS ,A5,I5,A5 , I5,A5,I5,A5,I5,A5 , 15 ,A5) 

777 FORMAT ( 2X , 7 I 5 , 1X,A8 , 8F10 . 5) 

199 FORMAT C5 14 ,8F7 .4) 

888 FORMAT (// ' NUMBER OF ELEMENTS IN THIS REGION =' ,2110/) 
400 CONTINUE 
C 


TRASFER 


PROGRAM 
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DIMENSION X(482) ,Y(482) ,T(482) ,BETA(376) ,PROB(3,7) ,ALFA(3,3) , 

1 BI (150) ,TA(482) ,TG(482) 

INTEGER NM , I , J , K , II , IE(376,5) , III , IT,MID(376) ,1G(482) , IS(482) 
CHARACTER*3 TIPE(376) ,PT(150) 

DOUBLE PRECISION PROB , ALFA , BETA 
DATA PROB /770000. 0,97500. 0,1000.0, 

1 164000.0,132600.0,1500.0, 

1 1380000.0,1950000.0,3000.0. 

1 0.440,0.287,0.40, 

1 0.166,0.022,0.02, 

1 0.255,0.02,0.02, 

1 76000.0,14000.0,500.0/ 

DATA ALFA /0 . 00002705 , 0 . 00003490 , 0 . 0000005 , 

1 0 . 00003038 , 0 . 0000259 , 0 . 0000005 , 

1 0.00001457,0.0000021,0.0000005/ 

DATA BIRT/'BIRT'/.BIRL/'BIRL’/, FI/' FI'/, PAP /‘PAP’/ 

DO 3 I * 1,3 

WRITE(7 ,55) (PROB(I , J) , J* 1,7) 

WRITE (7 ,44) (ALFA(I,J), J *1,3) 

CONTINUE 

FORMAT(3F10 . 1 , 3F10 . 8 ,F10 . 1) 

FORMAT (3F 10. 8) 

PI = 3.1415927 
TI = 70.0 


DO 5 I = 1,389 
IF(I.LE.140)TIPE(I) ='BIR' 

IF((I.GE.14l) .AND. (I.LE. 151)) TIPE(I) *’FIR' 
IF( (I .GE. 152) .AND. (I.LE.175)) TIPE(I) * , PAP l 
IFCCI.GE. 176). AND. (I.LE. 319)) TIPE(I) ='FIR' 
IFCCI.GE. 320). AND. (I.LE. 343)) TIPE(I) ='PAP' 
IF((I.GE. 344). AND. (I.LE. 376)) TIPE(I) *'FIR' 
IFCCI.GE. 377). AND. (I.LE. 389)) TIPE(I) ='BIR' 
CONTINUE 
DO 51 KK = 1,389 

IF (KK . GE . 1 . AND . KK . LE . 1 1 ) BETA (KK)=- . 26 18 
IF(KK.GE. 12.AND.KK.LE. 17) BETA(KK)=- . 15707 
IF (KK . GE . 18 . AND . KK . LE . 29 ) BETA (KK)=0 . 0 
IF (KK . GE . 30 . AND . KK . LE . 33 ) BETA (KK)=. 38394 
IF (KK . GE . 34 . AND . KK . LE . 37 ) BETA(KK)=. 17452 
IF (KK . GE . 38 . AND . KK. LE . 41) BETA (KK)=1. 22164 
IF (KK . GE . 42 . AND . KK . LE . 45 ) BETA (KK)=. 69808 
IF (KK . GE . 46 . AND . KK . LE . 49 ) BETA (KK)=2 .4433 
IF ( KK . GE . 5 0 . AND . KK . LE . 5 3 ) BETA(KK)=1.91972 

IFCKK.GE. 54.AND.KK.LE.57) BETA (KK)=2. 9668 
IFCKK.GE. 58.AND.KK.LE.61) BETA (KK)=2 . 7574 
IF (KK . GE . 62 . AND . KK . LE . 7 3 ) BETA (KK)=3 . 14136 
IF(KK.GE. 74.AND.KK.LE.79) BETA (KK)=3. 29843 
IFCKK.GE. 80 . AND . KK . LE . 10 1 ) BETA (KK)=3 .40314 
IFCKK.GE. 1G2 , AND . KK . LE . 107) BETA (KK)=3 . 29843 
IF (KK . GE . 108 . AND . KK . LE . 1 15 ) BETA(KK )=3 . 141 36 
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IF (KK . GE . 1 16 . AND . KK . LE . 1 23 ) BETA (KK)=0 . 0 
IF(KK . GE . 124 . AND . KK . LE . 129 ) BETA (KK) =-. 15707 
IF(KK.GE. 130 .AND.KK.LE. 163) BETA(KK)=-.2618 
IF (KK . GE . 1 64 . AND . KK . LE . 1 75 ) BETA (KK) =- . 1 5 7 0 7 
IF(KK.GE. 176. AND.KK.LE. 215) BETA(KK)=0.0 
IF(KK.GE, 216 .AND.KK.LE. 219) BETA(KK)=. 38394 
IF(KK . GE . 220 . AND . KK . LE . 223 ) BETA (KK)= . 17452 
IF(KK . GE . 224 . AND . KK. LE . 227 ) BETA (KK) = . 38394 
IF(KK.GE. 228. AND.KK.LE. 231) BETA (KK)=. 17452 
IF(KK.GE. 232. AND.KK.LE. 235) BETA(KK)=1.2216 
IF(KK.GE. 236. AND.KK.LE. 239) BETA (KK)=. 69808 
IFfKK. GE. 240. AND.KK.LE. 243) BETA (KK)=1 .2216 
IF (KK . GE . 244 . AND . KK . LE . 247 ) BETA (KK)“ . 69808 
IF (KK . GE . 248 . AND . KK . LE . 25 1 ) BETA (KK) =2 . 4433 
IF(KK.GE. 252. AND.KK.LE. 255) BETA(KK)=1 . 91972 
IF (KK . GE . 256 . AND . KK . LE . 25 9 ) BETA (KK) =2 . 4433 
IF (KK.GE. 260. AND.KK.LE. 263) BETA (KK )=1 . 91972 
IF(KK . GE . 264 . AND . KK . LE . 267 ) BETA (KK)=2 . 9668 
IF (KK . GE . 268 . AND . KK . LE . 27 1 j BETA (KK)=2 . 7574 
IF(KK . GE . 272 . AND . KK . LE . 275 ) BETA (KK)=2 . 9668 
IF(KK.GE. 276. AND.KK.LE. 279) BETA (KK)=2 . 7574 
IF (KK . GE , 280 . AND . KK . LE . 319 ) BETA (KK)=3 .14136 
IF(KK . GE . 320 . AND . KK . LE . 33 1 ) BETA (KK)=3 . 29843 
IF (KK . GE . 332 . AND . KK . LE . 354) BETA (KK)=3 . 403 14 
IF(KK.GE. 355. AND.KK.LE. 376) BETA(KK)=0.0 
C IF(KK.GE. 377. AND.KK.LE. 389) BETA (KK)=1 . 57079 
51 CONTINUE 

READ(5,200)(IG(I),I=1,482) 

C WRITE(6,200) (IG(I) ,1=1,482) 

READ(5, 200) (IS(I), 1=1,482) 

C WRITE(6,200) (IS(I) ,1=1,482) 

200 F0RMAT(5(6X,I5,5X)) 

DO 10 L = 1,376 

READ (5 ,777) NM, I , J,K,N,X(I) ,Y(I) ,X(J) ,Y(J) ,X(K) ,Y(K) ,X(N) ,Y(N) 
C WRITE (6, 777) NM, I , J,K,N,X(I) ,Y(I) ,X(J) ,Y(J) ,X(K) ,Y(K) ,X(N) ,Y(N) 

IE(L, 1) = NM 
IE(L, 2) = I 

IE(L,3) = J 
IE(L,4) = K 
IE(L,5) = N 
10 CONTINUE 

777 FORMAT (5 14 , 8F7 .4) 

READ (5 ,666) (T(I), I = 1,482) 

C WRITE (6, 666) (T(l), I = 1,482) 

DO 12 I = 1,482 
12 TG(IG(I))=T(IS(I)) 

666 FORMAT (9F8. 4) 

MID(l) =1 
DO 50 I = 2,376 

IF((TIPE(I) .NE.TIPE(I-l)) .OR. (BETA(I) .NE.BETA(I-l)))GO TO 25 
MID ( I ) = MID(I-l) 

GO TO 50 

•25 MID ( I ) = MID(I-l) + 1 

50 CONTINUE 
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IT = 0 

DO 70 I = 1,376 
IF (MID (I) .GE. IT) GO TO 15 
IT = MID (I) 

BI(IT) = BETA (I) 

PT(IT) = TIPE(I) 

WRITE (6, 16)1, TIPE(I) , BETA(I) , IT 
WRITE (6, 66) IT,PT(IT) 

CONTINUE 

F0RMAT(I5 ,5X,A4,F10.4, 15) 

DO 30 I = 1,482 
TA(I) = 90.0 + TG(I) 

WRITE (7, 555) I,X(I),Y(I) ,TA(I) 

CONTINUE 

F0RMAT(I5 ,40X,2F10 .4 ,5X,F10 . 4) 

DO 80 1=1, IT 
WRITE (6 ,66) I, PT(I) 

BII = BI(I)*180.0/PI 
WRITE(7,11) I, BII 
IF(PTCI).EQ.'BIR') GO TO 7 
IF(PT(I).EQ. 'FIR') GO TO 17 
IF(PT(I).EQ. 'PAP') GO TO 27 
WRITE (7, 22) TI, (PR0B(1, J) , J=l,7) 

WRITE(7,33) (ALFA(1,J), J =1,3) 

GO TO 80 

WRITE (7 ,22) TI , (PR0B(2 , J) , J=l,7) 

WRITE(7,33) (ALFA(2,J), J =1,3) 

GO TO 80 

WRITE (7 ,22) TI, (PR0B(3, J) , J=l,7) 

WRITE(7,33) (AEFA(3,J), J =1,3) 

CONTINUE 
FORMAT ( 15 , 5X , A4) 

FORMAT(I5,25X,F10.4) 

FORMAT (4F10 . 1 , 3F10 . 8 ,F10 .1) 

F0RMAT(3F10 .8) 

TH=1.000 
N= 4 

TR = 70.00 
M =1 

DO 40 I =1,376 

WRITE (7, 444) (IE(I, J) , J=l,5) ,MID(I) ,TR, N, TH 
CONTINUE 

FORMAT (6 15 ,F10 . 4 , 10X , 15 , 5X , F10 . 4) 

STOP 

END 
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14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 


4.3519 

13.0728 

147.8842 

3-8200 

12.5000 

144.3767 

3 . 2448 

14.0923 

150.2360 

2.6075 

13.3675 

147.7226 

5.0275 

13.6888 

150.4952 

4.0394 

14.8344 

152.0367 

3.4319 

11.9703 

139.9296 

2.1273 

12.6598 

144.3721 

2.1441 

15.1056 

154.4340 

1.4000 

14.2300 

154.0150 

3.0588 

15.9725 

154.7096 

0.8266 

13.3456 

153.5614 

5 . 8469 

14.3478 

152.3099 

4.9911 

15.5936 

153.2852 

4.1441 

16.8306 

154.9098 

3.1875 

11.4838 

134.7187 

1.8044 

11.9694 

139.7209 

0.4237 

12.4525 

152.5827 

2.0344 

15.2103 

154.7712 

1.2800 

14.3200 

154.5201 

2.9625 

16.0887 

154.9232 

0.6994 

13.4178 

154.2866 

4.0644 

16.9553 

155.0355 

0.2925 

12.5037 

153.5707 

6.8100 

15.0500 

153.5261 

6.1000 

16.3700 

154.1107 

5.4000 

17.6800 

155.0395 

5.3400 

17.8100 

155.1139 

3.0869 

11.0403 

129.1960 

1.6386 

11.2961 

133.0223 

0.1916 

11.5506 

150.4510 

0.0594 

11.5778 

154.5046 

7.2091 

15.3441 

153.9100 

6.5538 

16.6739 

154.3554 

5.9059 

18.0200 

155.0804 

5.8544 

18.1544 

155.1307 

3.1300 

10.6400 

124.1252 

1.6300 

10.6400 

123.5315 

0.1300 

10.6400 

123.5111 

0.0000 

10.6400 

123.9204 

7 . 8038 

15.6287 

154.2687 

7.2100 

16.9706 

154.5677 

6.6213 

18.3550 

155.1068 

6.5775 

18.4925 

155.1460 

3.0397 

10.2559 

119.1247 

1.5559 

10.0748 

115.2346 

0.0772 

9.8925. 

97.1668 


A 45 


* 


48 

-0.0594 

9.8781 

92.0556 

49 

8.5941 

15.9040 

154.5623 

50 

8.0688 

17.2601 

154.7453 

51 

7.5459 

18.6850 

155 . 1289 

52 

7.5094 

188243 

155.1588 

53 

3.0913 

9.8112 

113.2342 

54 

1.6412 

9.4569 

108.1028 

55 

0.2013 

9.1000 

93.2964 

56 

0.0625 

9.0725 

91.9792 

57 

9.5800 

16.1700 

154.7605 

58 

9.1300 

17.5425 

154.8743 

59 

8.6800 

19.0100 

155.1474 

60 

8.6500 

19.1500 

155.1694 

61 

3.2847 

9.3059 

107.4809 

62 

1.8859 

8.7861 

102.8803 

63 

0.5022 

8.2625 

91.9885 

64 

0 . 3656 

8.2231 

91.1432 

65 

10.7616 

16.4266 

154.8534 

66 

10.3938 

17.8177 

154.9407 

67 

10.0234 

19.3300 

155 . 1605 

68 

9.9994 

19.4693 

155.1775 

69 

3.6200 

8.7400 

102.4331 

70 

2.2900 

8.0625 

98.9507 

71 

0.9800 

7.3800 

91.4340 

72 

0.8500 

7.3300 

90.8650 

73 

12.1388 

16.6737 

154.8177 

74 

11.8600 

18.0856 

124.9131 

75 

11.5763 

19*6450 

155.1579 

76 

11.5575 

19.7825 

155.1759 

77 

4.0972 

8.1134 

98.3440 

78 

2.8534 

7.2861 

95.8928 

79 

1.6347 

6.4525 

90.9115 

80 

1.5156 

6.3931 

90.5528 

81 

13.7116 

16.9116 

154.5817 

82 

13.5288 

18.3464 

154.7357 

83 

13.3384 

19.9550 

155.1338. 

84 

13.3244 

20.0894 

155 . 1635 

85 

4.7163 

7.4263 

95.3039 

86 

3.5763 

6.4569 

93.6887 

87 

2.4663 

5.4800 

90.5535 

88 

2.3625 

5.4125 

90.3444 

89 

15.4800 

17.1400 

153.8904 

90 

15.4000 

18.6000 

154.2581 

91 

15.3100 

20.2600 

155.0557 

92 

15.3000 

20.3900 

155.1067 

93 

5.4772 

6.6784 

93.2149 

94 

4.4584 

5 .5748 

92.1772 

95 

3.4747 

4.4625 

90.3037 

96 

3.3906 

4.3881 

90.1936 

97 

16.3600 

17.2000 

153.0317 

98 

16.6300 

17.7400 

152.9341 

99 

15.8000 

17.8000 

153.7650 

100 

15.8025 

18.9350 

154.2794 

101 

15 . 8000 

20.2600 

155.0328 


A46 1 




102 

15.8000 

20.3900 

155.0989 

103 

6.3800 

5.8700 

91.8505 

104 

5.5000 

4.6400 

91.1968 

105 

4.6600 

3.4000 

90.1492 

106 

4.6000 

3.3200 

90.0924 

107 

16.8200 

17.2300 

152.0987 

108 

16.9700 

18.0400 

152.8914 

109 

17.0800 

18.8300 

153.7508 

110 

16.8200 

18.2800 

153.2712 

111 

16.3800 

18.3600 

153.6377 

112 

16.3400 

19.2200 

154.3084 

113 

16.3000 

20.2600 

155.0112 

114 

16.3000 

20.3900 

155.0851 

115 

6 . 7350 

5.5134 

91.4067 

116 

5.9581 

4.2261 

90.8833 

117 

5.2137 

2.9275 

90.1065 

118 

5.1578 

2.8425 

90.0704 

119 

17.2200 

17.2200 

150.9112 

120 

17.2200 

18.3000 

153.1369 

121 

17.2200 

19.3400 

154.2458 

122 

16.7600 

19.0000 

154.0092 

123 

17.2200 

19.4600 

154.3388 

124 

16.7600 

19.5750 

154.4896 

125 

16.7600 

20.2600 

154.9978 

126 

16.7600 

20.3900 

155.0772 

127 

7.2900 

5.1663 

90.9979 

128 

6.6150 

3.8294 

90.6383 

129 

5.9650 

2.4800 

90.0829 

130 

5.9137 

2.3900 

90.0547 

131 

17.8400 

17.5400 

151.4541 

132 

17.9050 

18.4900 

153.2845 

133 

17.2200 

15.6538 

144.1907 

134 

17.8725 

15.7737 

144.6731 

135 

17.9700 

19.3600 

154.2331 

136 

17.9700 

19.5000 

154.3391 

137 

17.2200 

19 . 8400 

154.6544 

138 

17.9700 

19.8750 

154.6646 

139 

17.2200 

20.2600 

154.9879 

140 

17.2200 

20.390 Q 

155.0715 

141 

8.0450 

4.8284 

90.6712 

142 

7.4706 

3.4498 

90.4450 

143 

6.9137 

2.0575 

90.0612 

144 

6.8678 

1.9625 

90.0401 

145 

18.4700 

17.8600 

152.5973 

146 

18.4700 

18.6700 

153.6167 

147 

18.4700 

15.8938 

145 .1270 

148 

18.4700 

19.3600 

154.2830 

149 

17.2200 

14.1600 

138,0302 

150 

17.9050 

14.1600 

138.0316 

151 

18.4700 

14.1600 

138.0318 

152 

18.4700 

19.5400 

154.4113 

153 

18.4700 

19.9100 

154.7148 

154 

17.9700 

20.2600 

154.9807 

155 

18.4700 

20.2600. 

154.9871 


156 

17.9700 

20.3900 

155.0665 

157 

9.0000 

4.5000 

90.4415 

158 

8.5250 

3.0875 

90.3062 

159 

8.0600 

1.6600 

90.0439 

160 

8.0200 

1.5600 

90.0285 

161 

19.4700 

17.8700 

153.7898 

162 

19.4700 

18.6550 

154.0946 

163 

19.4700 

19.3700 

154.5283 

164 

19.4700 

19.5200 

154.6021 

165 

17.2200 

12.7388 

132.1555 

166 

17.9375 

12.6987 

132.0249 

167 

18.4700 

12.6587 

131.8744 

168 

19.4700 

19.8925 

154.6245 

169 

19.4700 

20.2600 

155.0414 

170 

18.4700 

20.3900 

155.0680 

171 

19.4700 

20.3900 

155.1040 

172 

10.1550 

4.1809 

90.3129 

173 

9.7781 

2.7423 

90.2227 

174 

9.4038 

1.2875 

90.0314 

175 

9.3703 

1.1825 

90.0198 

176 

20.4700 

17,8800 

154.1075 

177 

20.4700 

18.6400 

154.3112 

178 

20.4700 

19.3800 

154.7841 

179 

20.4700 

19.5100 

154.8236 

180 

20.4700 

19.8800 

154.9594 

181 

17.2200 

11.3900 

126.2683 

182 

17.9700 

11.3900 

126.7595 

183 

18.4700 

11.3900 

126.8414 

184 

20.4700 

20.2600 

155.0976 

185 

20.4700 

20.3900 

155.1358 

186 

11.5100 

3.8713 

90.2962 

187 

11.2300 

2.4144 

90.2166 

188 

10.9450 

0.9400 

90.0313 

189 

10.9187 

0.8300 

90.0196 

190 

21.7062 

19.3500 

155.1174 

191 

21.7012 

19.4887 

155.1178. 

192 

21.7181 

19.8631 

155 . 1418 

193 

21.7325 

20.2350 

155 . 1785 

194 

17.2200 

10.6400 

124.3878 

195 

17.9700 

10.6400 

124.4096 

196 

16.4700 

10.6400 

124.3720 

197 

16.4700 

11.3900 

125.3677 

198 

18.4700 

10.6400 

124.4140 

199 

21.7387 

20.3550 

155.1893 

200 

13.0650 

3.5709 

90.4361 

201 

12.8806 

2.1036 

90.3264 

202 

12.6838 

0.6175 

90.0457 

203 

12.6653 

0.5025 

90.0261 

204 

22.9400 

19.3000 

155.1814 

205 

22 . 9400 

19.4400 

155.1815 

206 

22.9675 

19.8175 

155.1866 

207 

22.9900 

20.1800 

155. 1946 

208 

23.0000 

20.3000 

155.1970 

209 

17.2200 

9.8500. 

122.3997 
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210 

17.9700 

9.8500 

121.9564 

211 

16.4700 

9.8500 

123.2855 

212 

18.4700 

9.8500 

121.8808 

213 

14.8200 

3.2800 

90.9077 

214 

14.7300 

1.8100 

90.6688 

215 

14.6200 

0.3200 

90.1157 

216 

14,6100 

0.2000 

90.0768 

217 

24.1712 

19.2300 

155.1958 

218 

24.1862 

19.3637 

155.1959 

219 

24.2181 

19.7431 

155.1970 

220 

24 . 2425 

20.0950 

155.1988 

221 

24.2537 

20.2250 

155.1994 

222 

17.2200 

8,1750 

115.1947 

223 

17.9400 

8.2131 

115.3302 

224 

18.4700 

8.2537 

115.4882 

225 

15.2600 

2.7000 

91.1052 

226 

15.2650 

1.4925 

90.7231 

227 

16.0200 

3.2000 

91.7380 

228 

16.3000 

2.6500 

91.8801 

229 

15.2600 

0.2800 

90.1457 

230 

15 . 2400 

0.1500 

90.0878 

231 

25.4000 

19.1400 

155.1987 

232 

25.4400 

19.2600 

155.1988 

233 

25.4700 

19 . 6400 

155.1991 

234 

25.4900 

19.9800 

155.1997 

235 

25.5000 

20.1300 

155.1999 

236 

17.2200 

6.5000 

108.2843 

237 

17.9100 

6.4975 

108.2734 

238 

18.4700 

6.5000 

108.2838 

239 

15.9800 

2.1000 

91.3425 

240 

15.9400 

1.1600 

90.7478 

241 

16.6000 

2.1200 

91.7607 

242 

15.9000 

0.2300 

90.1710 

243 

16.6000 

3.1600 

92.6067 

244 

16.7350 

2.3650 

92.1014 

245 

16.8800 

1.5800 

91.3814 

246 

15.8800 

0.1000 

90; 1029 

247 

27.4722 

18.9144 

155.1998 

248 

27.5022 

19.0455 

155.1998 

249 

27.5397 

19.4275 

155.1999 

250 

27.5722 

19 . 7844 

155.2000 

251 

27.5855 

19.8989 

155.2000 

252 

17.2200 

4.8250 

101.3809 

253 

17.8800 

4.7031 

100.8879 

254 

18.4700 

4.5887 

100.4528 

255 

16.6000 

1.5100 

91.2262 

256 

16.5800 

0.8375 

90.6678 

257 

16.5600 

0.2000 

90.1982 

258 

16.5500 

0.0700 

90.1203 

259 

17.2200 

3.1500 

94.2286 

260 

17.2200 

2.1000 

92.0724 

261 

17.2200 

1.0500 

90.9383 

262 

17.2200 

0.9300 

90.8391 

263 

29.5188 

18.6544. 

155.2000 


Ps 


264 

29.5422 

18.7922 

155.2000 

265 

29.5856 

19.1733 

155.2000 

266 

29.6289 

19.5444 

155.2000 

267 

29.6455 

19.6355 

155.2000 

268 

17.8500 

2.8300 

93.6769 

269 

18.4700 

2.5200 

92.5586 

270 

17.2200 

0.5200 

90..4926 

271 

17.2200 

0.1700 

90.2098 

272 

17.2200 

0.0500 

90.1317 

273 

17.9100 

1.9350 

91.9496 

274 

17.9700 

1.0300 

90.9540 

275 

17.9700 

0.9100 

90.8619 

276 

17.9700 

0.5200 

90.5194 

277 

31.5400 

18.3600 

155.2000 

278 

31.5600 

18.5000 

155.2000 

279 

31.6075 

18.8775 

155.2000 

280 

31 . 6600 

19.2600 

155.2000 

281 

31.6800 

19.3400 

155.2000 

282 

18.4700 

1.7800 

91.6346 

283 

19.4700 

2.5400 

91.3945 

284 

19.4700 

1.7900 

91.1279 

285 

17.9700 

0.1600 

90.2209 

286 

17.9700 

0.0300 

90.1343 

287 

18.4700 

1.0200 

90.8996 

288 

18.4700 

0.9000 

90.8143 

289 

18.4700 

0.5200 

90.5006 

290 

18.4700 

0 . 1400 

90.2044 

291 

33.5355 

18.0311 

155.2000 

292 

33.5555 

18.1689 

155.2000 

293 

33.6055 

18.5400 

155.2000 

294 

33.6655 

18.9311 

155.2000 

295 

33.6888 

19.0122 

155.2000 

296 

19.4700 

1.0300 

90.6696 

297 

20.4700 

2.5600 

91.1005 

298 

20.4700 

1.8000 

90.9125 

299 

20.4700 

1.0400 

90.4442. 

300 

18.4700 

0 . 0200 

90.1300 

301 

19.4700 

0.9000 

90.6054 

302 

19.4350 

0.5150 

90.3788 

303 

19.4100 

0.1400 

90.1577 

304 

19.4100 

0.0200 

90.0994 

305 

35.5055 

17.6678 

155.2000 

306 

35.5289 

17.7989 

155.2000 

307 

35.5797 

18.1608 

155.2000 

308 

35 . 6455 

18.5578 

155.2000 

309 

35.6722 

18 . 6522 

155.2000 

310 

20.4700 

0.9100 

90.3963 

311 

21.2187 

1.0387 

90.1452 

312 

21.2175 

0.9162 

90.1430 

313 

20.4700 

0.5200 

90.2397 

314 

20.4900 

0.1500 

90.0937 

315 

20.4700 

0.0400 

90.0601 

316 

37.4500 

17.2700 

155.2000 

317 

37.4800 

17.3900. 

155.2000 


A 


316 

37.5299 

17.7400 

155.2000 

319 

37.6000 

18.1400 

155.2000 

320 

37.6300 

18.2600 

155.2000 

321 

21.2100 

0.5300 

90.1030 

322 

22.2900 

1.0800 

90.0357 

323 

22.2900 

0.9600 

90.0355 

324 

22.2850 

0.5800 

90.0275 

325 

21.2175 

0.1538 

90.0439 

326 

21.2100 

0.0375 

90.0284 

327 

37 . 8422 

17.2166 

155.2000 

328 

37.8555 

17 . 3455 

155.2000 

329 

37.8680 

17.7036 

155.2000 

330 

37.8988 

18.1066 

155.2000 

331 

37.9122 

18.2377 

155.2000 

332 

22.2900 

0.2000 

90.0110 

333 

23.6837 

1.1637 

90.0063 

334 

23.6875 

1.0413 

90.0062 

335 

23.6950 

0.6700 

90.0048 

336 

23.7075 

0.2887 

90.0019 

337 

22.2900 

0.0800 

90.0062 

338 

38.8589 

16.9933 

155.2000 

339 

38.8622 

17.1289 

155.2000 

340 

38.8555 

17.4911 

155.2000 

341 

38.8655 

17.8933 

155.2000 

342 

38.8689 

18.0311 

155.2000 

343 

23.7100 

0.1675 

90.0010 

344 

25.4000 

1.2900 

90.0013 

345 

25.4100 

1.1600 

90.0013 

346 

25.4400 

0.8000 

90.0009 

347 

25.4700 

0.4200 

90.0002 

348 

25.4700 

C.3000 

90.0001 

349 

40.5000 

16.6000 

155.2000 

350 

40.5000 

16.7400 

155.2000 

351 

40.4925 

17.1025 

155.2000 

352 

40.5000 

17.5000 

155.2000 

353 

40.5000 

17.6400 

155.2000 

354 

27.3933 

1.5367 

90.0002 

355 

27.4111 

1.4111 

90.0002 

356 

27.4577 

1.0478 

90.0001 

357 

27.5011 

0.6711 

90.0000 

358 

27.5155 

0.5456 

90.0000 

359 

42.7655 

16.0367 

155.2000 

360 

42.7689 

16.1789 

155.2000 

361 

42.7788 

16.5378 

155.2000 

362 

42.8022 

16.9267 

155.2000 

363 

42.8055 

17.0644 

155.2000 

364 

29.3733 

1.8133 

90.0000 

365 

29.3978 

1.6911 

90.0000 

366 

29.4577 

1.3244 

90,0000 

367 

29.5111 

0.9511 

90.0000 

368 

29.5355 

0.8222 

90.0000 

369 

45 .6555 

15 . 3033 

155.2000 

370 

45.6689 

15.4456 

155.2000 

371 

45.7147 

15.7970. 

155.2000 


372 

45.7722 

16.1733 

155 . 2000 

373 

45 . 7855 

16.3044 

155.2000 

374 

31.3400 

2.1200 

90.0000 

375 

31.3700 

2.0000 

90.0000 

376 

31,4400 

1.6300 

90.0000 

377 

31.5000 

1.2600 

90.0000 

378 

31.5300 

1 . 1300 

90.0000 

379 

49.1700 

14.4000 

155.2000 

380 

49.2000 

14.5400 

155 . 2000 

381 

49.2999 

14.8800 

155.2000 

382 

49.4100 

15.2400 

155.2000 

383 

49.4400 

15.3600 

155.2000 

384 

33.2933 

2.4567 

90.0000 

385 

33.3277 

2.3378 

90.0000 

386 

33.4044 

1.9644 

90.0000 

387 

33.4677 

1.5978 

90.0000 

388 

33.4989 

1.4689 

90.0000 
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212 

210 

44 

70.0000 

4 

1.0000 

357 

236 

237 

223 

222 

44 

70.0000 

4 

1.0000 

358 

237 

238 

224 

223 

44 

70.0000 

4 

1.0000 

359 

252 

253 

237 

236 

44 

70.0000 

4 

1.0000 

360 

253 

254 

238 

237 

44 

70.0000 

4 

1.0000 

361 

259 

268 

253 

252 

44 

70.0000 

4 

1.0000 

362 

268 

269 

254 

253 

44 

70.0000 

4 

1.0000 

363 

194 

195 

182 

181 

44 

70.0000 

4 

1.0000 

364 

195 

198 

183 

182 

44 

70.0000 

4 

1.0000 

365 

209 

210 

195 

194 

44 

70.0000 

4 

1.0000 

36 6 

210 

212 

198 

195 

44 

70.0000 

4 

1.0000 

367 

196 

194 

181 

197 

44 

70.0000 

4 

1.0000 

368 

211 

209 

194 

196 

44 

70.0000 

4 

1.0000 

369 

133 

134 

131 

119 

44 

70.0000 

4 

1.0000 

370 

134 

147 

145 

131 

44 

70.0000 

4 

1.0000 

371 

149 

150 

134 

133 

44 

70.0000 

4 

1.0000 

372 

150 

151 

147 

134 

44 ' 

70.0000 

4 

1.0000 

373 

165 

166 

150 

149 

44 

70.0000 

4 

1.0000 

374 

166 

167 

151 

150 

44 

70.0000 

4 

1.0000 

375 

181 

182 

166 

165 

44 

70.0000 

4 

1.0000 

376 

182 

183 

167 

166 

44 

70.0000 

4 

1.0000 

377 

434 

485 

484 

435 

45 

70 . 0000 

4 

1.0000 

378 

433 

486 

485 

434 

45 

70.0000 

4 

1.0000 

379 

430 

487 

486 

433 

45 

70.0000 

4 

1.0000 

380 

429 

488 

487 

430 

45 

70.0000 

4 

1.0000 

381 

432 

489 

488 

429 

45 

70.0000 

4 

1.0000 

382 

483 

490 

489 

432 

45 

70.0000 

4 

1.0000 

383 

478 

491 

490 

483 

45 

70.0000 

4 

1 ,0000 


At 63 


384 

476 

492 

491 

478 

45 

70.0000 

4 

385 

477 

493 

492 

476 

45 

70.0000 

4 

386 

480 

494 

493 

477 

45 

70.0000 

4 

387 

481 

495 

494 

480 

45 

70.0000 

4 

388 

482 

496 

495 

481 

45 

70.0000 

4 


/* 

// 


1.0000 

1.0000 

1.0000 

1.0000 

1.0000 


1.0 
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RADIAL-TANGENTIAL STRESS PRGRAM 

DIMENSION S2ETA(90) .SEATAC90) ,SZT(90) ,BETA(90) 

REAL Sll, S22, S12, AA, BB 
INTEGER NN, NEL 
DO 20 KK = 1,90 

IFCKK.GE. 1. AND.KK.LE. 11) BETACKK)=-.26l8 
IFCKK.GE. 12.AND.KK.LE. 17) BETA (KK)=-. 1S707 
IF(KK. GE . IS . AND . KK . LE . 29 ) BETA(KK)=0 . 0 
IF (KK . GE . 30 . AND . KK . LE . 33 ) BETA (KK)= . 38394 
IF (KK . GE . 34 . AND , KK . LE . 37 ) BETA (KK)=. 17452 
IF (KK . GE . 38 . AND . KK . LE . 4 1 ) BETA (KK)= 1.22164 
IF (KK . GE . 42 . AND . KK. LE . 45 ) BETA (KK)= .69808 
IF (KK . GE . 46 . AND . KK . LE . 49 ) BETA (KK)=2 .4433 
IF(KK.GE. 50 .AND.KK.LE.53) BETA(KK)=1 .91972 
IFCKK.GE. 54.AND.KK.LE.57) BETA (KK)=2. 9668 
IFCKK.GE. 58 .AND.KK.LE. 61) BETA(KK)=2 . 7574 
IFCKK.GE. 62. AND.KK.LE. 73) BETA (KK)=3. 14136 
IFCKK.GE. 74. AND.KK.LE. 79) BETA CKK)=3 .29843 
IFCKK.GE. 80. AND.KK.LE. 90) BETA CKK) =3 . 40314 
20 CONTINUE 

READ C5,11)NN 
11 FORMAT (1 5) 

WRITE C6, 15) 

15 FORMAT C / / , 1 OX , 3HNEL , 9X , 6HS - ZETA , 9X , 6HS - E ATA , 9 X , 10HTUE - ZET-ET,//) 

DO 10 1=1, NN 

RE AD C5, 22) NEL, Sll, S22, S12 

AA=CSll+S22)/2 

BB= (Sll-S22)/2 

SZETA(NEL)=AA+BB*C0S(2*BETA(NEL))+S12*SIN(2*BETA(NEL)) 
SEATA(NEL)=AA-BB*C0S(2*BETACNEL)) -S12*SIN(2*BETA(NEL) ) 
SZTCNEL) =-BB*SINC2*BETACNEL) )+S12*C0S(2*BETACNEL)) 

WRITE C6 , 33)NEL, SZETA(NEL) , SEATACNEL) , SZTCNEL) 

10 CONTINUE 

22 FORMAT C 15 , 3E15 . 5) 

33 FORMAT (5X , 15 , 5X , 3F15 . 3 , / ) 

STOP 

END 


A <>5 


TEMPERATURE IN DEGREE F 


NODE NO. TEMP. NODE NO. TEMP. 


1 

64.771 

2 

64.015 

5 

64.923 

6 

63.561 

9 

64.710 

10 

60 . 236 

13 

62.583 

14 

49 . 721 

17 

54.377 

18 

64.910 

21 

63.571 

22 

65.114 

25 

39.196 

26 

65.039 

29 

60.495 

30 

64.505 

33 

33.511 

34 

33.532 

37 

62.310 

38 

63.526 

41 

65.107 

42 

65 . 146 

45 

29 . 125 

46 

64.568 

49 

65 . 129 

50 

65.159 

53 

23.234 

54 

64.745 

57 

65 . 147 

58 

65.169 

61 

12.880 

62 

17.481 

65 

65.161 

66 

65.177 

69 

8.951 

70 

12.433 

73 

65 . 158 

74 

65.176 

77 

5.893 

78 

8.344 

81 

65 . 134 

82 

65.163 

85 

3.689 

86 

5.304 

89 

65.056 

90 

65.107 

93 

0.194 

94 

2.177 

97 

65.033 

98 

65.099 

101 

0.149 

102 

1.197 

105 

63.032 

106 

65.011 

109 

63.638 

110 

62.934 

113 

1.407 

114 

0.070 

117 

65.077 

118 

64.490 

121 

62.891 

122 

0.083 

125 

0.055 

126 

60.911 

129 

64.654 

130 

64.339 

133 

63.137 

134 

0.061 

137 

0.040 

138 

61.454 

141 

64.981 

142 

65.067 

145 

64.233 

146 

63 . 285 

149 

0.441 

150 

0.029 

153 

48.030 

154 

48.032 

157 

65 . 068 

158 

64.715 

161 

63.617 

162 

0.031 

165 

0.020 

166 

63.790 

169 

42.025 

170 

41.874 

173 

64.824 

174 

64.602 

177 

0.217 

178 

0.296 

181 

64.311 

182 

64.784 

185 

36.841 

186 

65.098 


NODE NO. TEMP. NODE NO. TEMP. 


3 

64.434 

4 

64.520 

7 

54.372 

8 

57.723 

11 

64.287 

12 

65.035 

15 

44.719 

16 

49.930 

19 

62.037 

20 

57.884 

23 

60.451 

24 

43.022 

27 

63.285 

28 

64.111 

31 

65.080 

32 

65.131 

35 

34.125 

36 

64.355 

39 

63.910 

40 

33.920 

43 

7.167 

44 

25.235 

47 

64.269 

48 

2.056 

51 

3.296 

52 

18.103 

55 

64.562 

56 

1.979 

59 

1.989 

60 

1.143 

63 

64.874 

64 

64.760 

67 

1.434 

68 

0.865 

71 

64.941 

72 

64.853 

75 

0.911 

76 

0.553 

79 

64.913 

80 

64.818 

83 

0.554 

84 

0.344 

87 

64.736 

88 

64.582 

91 

64.258 

92 

0.304 

95 

3.215 

96 

63.890 

99 

64.279 

100 

63.765 

103 

0.092 

104 

1.851 

107 

65.085 

108 

64.308 

111 

0.107 

112 

0.883 

115 

62.099 

116 

64.998 

119 

64.009 

120 

63.271 

123 

0.638 

124 

0.998 

127 

64.988 

128 

65.071 

131 

63.751 

132 

64.246 

135 

0.445 

136 

0.671 

139 

54.191 

140 

54.673 

143 

64.665 

144 

64.339 

147 

0.044 

148 

0.306 

151 

62.597 

152 

55.127 

155 

48.032 

156 

64.987 

159 

64.411 

160 

64.283 

163 

0.223 

164 

0.313 

167 

64.095 

168 

42.156 

171 

65.041 

172 

65 . 104 

175 

64.528 

176 

0.031 

179 

0.020 

180 

64.107 

183 

36.268 

184 

36.760 

187 

65.136 

188 

.64.959 


189 

64.824 

190 

0.046 

191 

0.326 

192 

0.436 

193 

0.026 

194 

65.117 

195 

65.118 

196 

34.388 

197 

34.372 

198 

35.368 

199 

34.410 

200 

34.414 

201 

65.179 

202 

65.189 

203 

65.142 

204 

0.116 

205 

0.669 

206 

0.908 

207 

0.077 

208 

65.181 

209 

65.182 

210 

65.187 

211 

32.400 

212 

33.285 

213 

31.956 

214 

31.881 

215 

65.195 

216 

65.197 

217 

0.146 

218 

0.723 

219 

1.105 

220 

1.738 

221 

0.088 

222 

65.196 

223 

65.196 

224 

65.197 

225 

65 . 199 

226 

25.195 

227 

25.330 

228 

25.488 

229 

65.199 

230 

0.171 

231 

0.103 

232 

0.748 

233 

1.343 

234 

1.880 

235 

2.607 

236 

65 . 199 

237 

65.199 

238 

65.199 

239 

65.200 

240 

65,200 

241 

18.284 

242 

18.273 

243 

18.284 

244 

0.198 

245 

0.120 

246 

0 . 668 

247 

1.226 

248 

1.761 

249 

2.101 

250 

4.229 

251 

65.200 

252 

65 . 200 

253 

65.200 

254 

65.200 

255 

65.200 

256 

11.381 

257 

10.888 

258 

10.453 

259 

0.210 

260 

0.493 

261 

0.132 

262 

0.839 

263 

1.381 

264 

2.072 

265 

3.677 

266 

65.200 

267 

65.200 

268 

65.200 

269 

65.200 

270 

65.200 

271 

2.559 

272 

0.221 

273 

0.519 

274 

0.862 

275 

0.134 

276 

0.954 

277 

0.938 

278 

1.950 

279 

65.200 

280 

65.200 

281 

65.200 

282 

65.200 

283 

65.200 

284 

1.635 

285 

1.394 

286 

0.204 

287 

0.501 

288 

0.814 

289 

0.900 

290 

0.130 

291 

65.200 

292 

65.200 

293 

65.200 

294 

65.200 

295 

65.200 

296 

1.128 

297 

1.100 

298 

0.158 

299 

0.379 

300 

0.670 

301 

0.605 

302 

0.099 

303 

65.200 

304 

65 . 200 

305 

65.200 

306 

65.200 

307 

65.200 

308 

0.913 

309 

0.094 

310 

0.060 

311 

0.240 

312 

0.444 

313 

0.397 

314 

65.200 

315 

65.200 

316 

65.200 

317 

65.200 

318 

65.200 

319 

0.044 

320 

0.028 

321 

0.103 

322 

0.143 

323 

0.145 

324 

65.200 

325 

65.200 

326 

65.200 

327 

65.200 

328 

65 . 200 

329 

0.011 

330 

0.028 

331 

0.006 

332 

0.036 

333 

0.035 

334 

65.200 

335 

65.200 

336 

65 . 200 

337 

65.200 

338 

65.200 

339 

0.002 

340 

0.005 

341 

0.001 

342 

0.006 

343 

0.006 

344 

65.200 

345 

65.200 

346 

65.200 

347 

65.200 

348 

65.200 

349 

0.000 

350 

0.001 

351 

0.001 

352 

0.000 

353 

0.001 

354 

65.200 

355 

65.200 

356 

65.200 

357 

65.200 

358 

65.200 

359 

0.000 

360 

0.000 

361 

0.000 

362 

0.000 

363 

- 0.000 

364 

65.200 

365 

65.200 

366 

65.200 

367 

65.200 

368 

65.200 

369 

0.000 

370 

0.000 

371 

0.000 

372 

0.000 

373 

0.000 

374 

65.200 

375 

65 . 200 

376 

65.200 

377 

65 . 200 

378 

65.200 

379 

- 0.000 

380 

0.000 

381 

0.000 

382 

0.000 

383 

- 0.000 

384 

65.200 

385 

65.200 

386 

65.200 

387 

65.200 

388 

65 . 200 

389 

0.000 

390 

0.000 

391 

0.000 

392 

0.000 

393 

0.000 

394 

65.200 

395 

65.200 

396 

65.200 

397 

65 . 200 

398 

65.200 

399 

- 0.000 

400 

0.000 

401 

0.000 

402 

0.000 

403 

- 0.000 

404 

.65.200 
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405 

65 . 200 

406 

65.200 

407 

65 . 200 

408 

65.200 

409 

0.000 

410 

0.000 

411 

0.000 

412 

0.000 

413 

0.000 

414 

65.200 

415 

65.200 

416 

65.200 

417 

65.200 

418 

65 , 200 

419 

- 0.000 

420 

0.000 

421 

0.000 

422 

0.000 

423 

- 0.000 

424 

65.200 

425 

65.200 

426 

65.200 

427 

65.200 

428 

65 . 200 

429 

65.200 

430 

0.000 

431 

0.000 

432 

0.000 

433 

0.000 

434 

0.000 

435 

65.200 

436 

- 0.000 

437 

0.000 

438 

0.000 

439 

0.000 

440 

- 0.000 

441 

0.000 

442 

0.000 

443 

0.000 

444 

0.000 

445 

0.000 

446 

- 0.000 

447 

- 0.000 

448 

0.000 

449 

0.000 

450 

- 0.000 

451 

0.000 

452 

0.000 

453 

0.000 

454 

0.000 

455 

0.000 

456 

0.000 

457 

- 0.000 

458 

0.000 

459 

0.000 

460 

0.000 

461 

0.000 

462 

0.000 

463 

0.000 

464 

0.000 

465 

0.000 

466 

0.000 

467 

- 0.000 

468 

0.000 

469 

0.000 

470 

0.000 

471 

0.000 

472 

0 . 000 

473 

0.000 

474 

0.000 

475 

0.000 

476 

0.000 

477 

0.000 

478 

0.000 

479 

0.000 

480 

0.000 

481 

0.000 

482 

0.000 
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NEL 

°Z 

S-ZETA 

1 

72.753 

2 

128.730 

3 

150.471 

4 

134.717 

S 

103.118 

6 

42.985 

7 

43.985 

8 

40.778 

9 

36.630 

10 

34.134 

n 

30.389 

12 

58.134 

13 

71.830 

14 

84.303 

15 

85.364 

16 

74.380 

17 

57.161 

18 

182.210 

19 

193.090 

20 

184.950 

21 

138.810 

22 

141.180 

23 

159.040 

24 

159.540 

25 

143.690 

26 

181.040 

27 

175.540 

28 

206.120 

29 

196.080 

30 

219.495 

31 

220.655 

32 

219.076 

33 

218.126 

34 

220.024 

35 

219.586 

36 

219.361 

37 

216.837 

38 

212.754 

39 

281.875 

40 

228.291 

41 

198.368 

42 

196.845 

43 

205.437 

44 

211.376 

45 

218.019 

46 

54.648 

47 

60.289 





S-EATA 

TUE-ZET-ET 


-25.031 

0.097 

21.610 

-10.909 

17.487 

-3.885 

7.048 

-3.030 

-6.019 

-6.335 

1.648 

1.088 

8.558 

0.592 

4.444 

0.450 

-2.692 

0.738 

-8.519 

-4.441 

-9.585 

6.986 

-24.751 

-12.546 

-7.608 

4.800 

-0.846 

1.810 

4.561 

-1.291 

5.382 

-0.979 

0.881 

-0.938 

6.823 

-4.466 

13.623 

-6.950 

17.615 

-4.893 

5.602 

2.921 

-6.364 

3.446 

-0.724 

8.298 

-1.576 

-4.006 

-5.850 

-3.097 

-1.726 

-2.301 

1.292 

-2.021 

-2.201 

-0.493 

-0.605 

-2.915 

2.378 

8.540 

3.084 

3.971 

0.208 

-4.697 

-1.217 

-8.853 

-0.417 

2.685 

1.129 

3.916 

0.641 

0.488 

-0.101 

-3.518 

62.153 

24.500 

95.675 

-25.342 

42.228 

-19,808 

4.434 

-18.915 

-7.460 

11.578 

1.011 

6.168 

0.407 

-3.798 

-1.705 

-8.118 

-2,145 

-4.063 

-0.573 

-6.099. 
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48 

64.369 

0.498 

- 9.703 

49 

70.544 

2.808 

- 12.586 

50 

71.693 

- 12.204 

- 2.760 

51 

44.375 

- 42.569 

- 23.700 

52 

- 14.116 

- 79.770 

- 49.946 

53 

21.449 

- 105.675 

- 27.785 

54 

- 0.499 

- 1.587 

- 1.159 

55 

- 9.181 

- 0.671 

- 2.958 

56 

18.219 

0.743 

- 0.578 

57 

25.410 

0.581 

0.050 

58 

34.978 

1.825 

- 0.734 

59 

40.194 

1.519 

- 1.884 

60 

39.774 

1.679 

- 5.518 

61 

32.715 

- 3.919 

- 8.640 

62 

- 9.067 

- 1.591 

1.703 

63 

- 6.920 

- 3.156 

1.557 

64 

- 13.227 

- 0.119 

0.270 

65 

- 11.307 

- 0.518 

1.700 

66 

- 9.902 

2.126 

- 1.099 

67 

- 11.378 

0.579 

- 2.753 

68 

- 33.096 

- 5.307 

3.648 

69 

- 15.536 

1.261 

2.706 

70 

- 36.337 

6.309 

- 2.919 

71 

- 40.329 

0.062 

- 1.569 

72 

- 44.515 

- 2.928 

- 0.573 

73 

- 47.240 

- 4.734 

- 0.251 

74 

- 55.141 

0.136 

- 3.032 

75 

- 55,361 

- 1.528 

- 1.849 

76 

- 56.541 

- 0.846 

- 0.927 

77 

- 58.733 

0.498 

- 0.432 

78 

- 62.359 

1.659 

0.193 

79 

- 67 . 046 

- 0.531 

5.309 

80 

- 32.720 

- 6.246 

- 3.719 

81 

- 38,555 

- 5.031 

0.494 

82 

- 44.718 

- 3.183 

0.305 

83 

- 50.290 

- 0.961 

0.029 

84 

- 54.344 

2.025 

0.268 

85 

- 56.160 

4.501 

1.446 

86 

1.913 

2.701 

- 3.104 

87 

1.889 

1.905 

- 0.896 

88 

1.010 

1.776 

0.904 

89 

- 5.577 

0.805 

3.054 

90 

1.411 

- 4.648 

- 8.905 
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